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Executive  Summary 


Explosions  in  mixtures  of  natural  gas  (NG)  and  air  have  been  of  intense  practical  concern  for  coal  mines  for 
many  years.  Potentially  explosive  mixtures  of  NG  and  air  can  accumulate  in  the  active  ventilated  areas  or 
in  unventilated  sealed  areas  of  these  mines.  If  an  ignition  source,  such  as  a  simple  spark,  ignites  the  NG- 
air  mixture  and  creates  a  flame,  the  initially  slow-moving  flame  can  become  turbulent,  accelerate  rapidly, 
develop  extremely  intense  pressure  waves,  and  potentially  generate  enormous  stress  on  coal  mine  seals. 
Recent  NIOSH  and  US  ACE  studies  of  an  explosion  at  the  Sago  mine  suggested  that  high  pressures  and 
dcflagration-to-detonation  transition  (DDT)  can  occur  in  NG-air  mixtures  in  coal  mines,  and  that  pressure 
criteria  for  seal  design  may  be  inadequate. 

In  this  project  sponsored  by  NIOSH-OMSHR,  research  efforts  by  NRL  focused  on  a  critical,  overriding 
question  for  coal-mine  safety:  Given  a  large  enough  volume  of  a  flammable  mixture  of  NG  and  air,  such  as 
may  exist  in  a  coal  mine,  can  a  weak  spark  ignition  develop  into  a  detonation? 

To  address  this  question,  we  developed  computational  models  based  on  the  reactive  Navier-Stokes 
equations  (RNSE)  that  can  be  used  to  compute  combustion  waves  ranging  from  low-speed  laminar  flames, 
to  moderate-speed  turbulent  flames,  to  high-speed  detonations.  The  combination  of  efficient  numerical 
algorithms,  adaptive  mesh  refinement,  and  one-step  reaction  kinetics  used  in  these  models  allows  a  wide 
range  of  scales  to  be  resolved  in  multidimensional  computations.  The  model  parameters  were  calibrated 
to  reproduce  known  properties  of  laminar  flames  and  detonations  in  various  methane-air  mixtures,  and  the 
RNSE  simulations  were  validated  using  experimental  data  on  DDT  in  obstructed  channels.  The  numerical 
models  were  used  to  study  four  specific  problems: 

1.  Flame  acceleration  and  DDT  in  channels  containing  a  stoichiometric  methane-air  mixture 
and  obstacles. 

We  have  shown  that  it  is  possible  to  compute  DDT  in  nrethane-air  mixtures  from  first  principles  by  solv¬ 
ing  the  RNSE,  given  enough  numerical  resolution  of  critical  regions  in  the  system  and  a  chemical  model 
that  is  valid  for  the  combustion  regimes  ranging  from  autoignition  to  laminar-flame  propagation.  Flames 
in  obstructed  channels  accelerate  to  supersonic  velocities  and  produce  shocks.  Detonations  appear  when 
shocks  interacting  with  obstacles  become  strong  enough  to  ignite  the  unburned  material.  Simulations 
performed  for  a  number  of  experimental  configurations  show  good  agreement  between  numerical  and 
experimental  results,  and  thus  they  can  be  used  as  a  predictive  tool. 

2.  Flame  acceleration  and  DDT  in  fuel-rich  and  fuel-lean  mixtures  in  channels  with  obstacles. 

Simulations  performed  for  rich  and  lean  mixtures  showed  that  the  physical  events  leading  to  DDT  in 
obstructed  channels  are  essentially  the  same  as  those  in  stoichiometric  mixtures.  The  results  show  a 
reasonable  agreement  with  a  limited  number  of  available  experiments. 

3.  Effects  of  spatially  varying  fuel  concentrations  on  detonation  propagation. 

A  single-step,  two-component  reaction  model  for  nonuniform  methane-air  mixtures  was  developed  and 
used  in  simulations  of  a  detonation  propagating  through  the  concentration  gradient  in  a  direction  normal 
to  this  gradient.  The  computed  reaction-zone  structure  involved  two  distinct  burning  regimes.  The 
leading  edge  of  the  reaction  zone  is  an  unstable  detonation  wave  with  a  variable  cell  size  that  depends 
on  a  local  mixture  composition.  Behind  the  detonation,  the  excess  fuel  and  oxidizer  that  are  not 
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consumed  by  detonations  in  rich  and  lean  mixtures,  respectively,  mix  and  burn  in  a  turbulent  diffusion 
flame.  A  series  of  independent  calculations  of  detonations  and  flames  have  been  performed  to  test  and 
validate  the  single-step,  two-component  reaction  model.  We  plan  to  use  the  same  reaction  model  to 
study  flame  acceleration  and  DDT  in  mixtures  with  a  nonunifornr  composition. 

4.  Effects  of  stochasticity  on  flame  acceleration  and  DDT. 

Experimental  observations  show  that  DDT  is  a  stochastic  process  in  the  sense  that  there  is  some  uncer¬ 
tainty  in  time  and  location  for  the  detonation  initiation.  In  simulations  of  DDT  in  obstructed  channels, 
the  onset  of  DDT  is  also  found  to  vary  stochastically,  and  its  location  is  only  predictable  with  an 
uncertainty  15-20%,  which  is  about  the  same  as  observed  in  experiments. 

The  numerical  studies  were  performed  by  NRL  in  conjunction  with  experimental  work  by  NIOSH.  The 
experiments  were  performed  in  the  Gas  Explosion  Test  Facility  (GETF),  a  new  1.05  m  by  72  m  detonation 
tube  at  the  NIOSH  Lake  Lynn  Laboratory.  Experiments  performed  in  GETF  have  now  extended  the  known 
detonability  limits  for  methane-air  mixtures  to  the  include  mixtures  of  5.3%  to  15.5%  methane  in  air.  They 
have  also  shown  that  DDT  can  occur  in  the  1.05  m  diameter  tube.  The  NRL  team  provided  expert  support 
for  the  design  and  implementation  of  this  facility. 

Both  computational  and  experimental  studies  have  made  substantial  progress  toward  answering  the 
fundamental  question:  Is  DDT  and  detonation  possible  in  coal  mine  tunnels?  To  provide  a  definitive, 
scientifically  sound  answer,  we  still  need  to  solve  three  critical  problems: 

1.  We  must  prove  that  numerical  models  calibrated  on  small-scale  experiments  can  be  extrapolated  to  the 
larger  geometry  of  typical  coal  mine  tunnels.  Additional  validation  tests  with  the  1.05  m  tube  at  NIOSH 
Lake  Lynn  Laboratory  coupled  with  a  few  large  scale  tests  at  Lake  Lynn  Experimental  Mine  should 
provide  the  critical  validation  data. 

2.  We  must  understand  and  predict  flame  behavior,  DDT  and  detonation  in  spatially  nonunifornr  methane- 
air  mixtures  that  are  likely  to  form  in  real-life  situations. 

3.  We  must  understand  the  mechanism  for  DDT  in  smooth  channels,  which  has  been  observed  in  ex¬ 
periments  but  has  not  been  computed  with  current  models,  and  remains  a  key  unresolved  problem  in 
combustion. 
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1.  Introduction 


Explosions  in  mixtures  of  natural  gas  (NG)  and  air  have  been  the  subject  of  intense  practical  concern  for 
many  years.  Such  explosions  occur  naturally  in  coal  mines  where  NG  accumulates,  and  accidentally  when 
vessels  or  pipelines  containing  NG  or  liquid  NG  are  ruptured.  Thus,  maintaining  safe  conditions  in  mines 
and  in  the  processes  of  extracting,  handling,  and  transporting  NG  is  of  major  importance  to  our  health  and 
economy.  Care  must  therefore  be  taken  to  define  and  avoid  conditions  which  can  lead  to  these  explosions. 
There  are  many  reviews  and  references  to  these  types  of  large-scale  explosions  (see,  for  example,  [1-4]),  and 
these  explicitly  show  the  explosive  dangers  of  NG  and  other  fuels. 

Large-scale,  gaseous  explosions,  either  natural  or  accidental,  lead  to  destructive  scenarios  for  which  we 
pay  heavily  in  lives  and  property.  Consider  some  of  the  recent  large-scale  events  that  we  have  read  and 
heard  about  in  the  media.  In  2006,  the  explosion  in  the  Sago  Mine  trapped  thirteen  miners  for  two  days, 
and  only  one  survived.  That  same  year,  there  was  an  explosion  in  a  fuel  storage  plant  in  Buncefield,  UK, 
which  was  heard  as  far  as  20  miles  away.  Because  this  explosion  occurred  late  at  night,  there  were  only  few 
people  in  the  immediate  vicinity  and,  miraculously,  no  loss  of  life.  The  cost  in  property  and  damages  in 
the  plant  and  surrounding  houses  exceeded  a  billion  dollars.  In  2010,  there  was  a  large  explosion,  heard  ten 
miles  away,  in  a  power  plant  being  constructed  in  Middletown,  Connecticut.  This  cost  lives,  injured  many, 
and  cost  millions  of  dollars  in  damages. 

NG  explosions,  are  extremely  complex  events.  In  coal  mines,  explosions  can  result  from  ignition  of 
NG  accumulated  in  unventilated,  sealed  areas.  In  time,  the  NG  mixes  with  the  ambient  air  and  forms  a 
flammable  mixture.  It  might  only  take  a  small  spark  to  ignite  NG  and  create  a  flame,  which  then  interacts 
with  its  environment,  accelerates,  generates  high  pressures,  and  may  eventually  lead  to  a  detonation.  Figure 
1.1  (which  is  discussed  more  thoroughly  in  Section  3  of  this  report)  shows  how  this  can  occur  in  a  methane-air 
mixture.  The  small  flame  naturally  becomes  turbulent  as  it  generates  pressure  waves  that  reflect  from  walls  of 
the  confining  chamber.  Eventually,  shock  waves  arise  as  these  waves  coalesce,  and  these  continue  to  interact 
with  the  walls  and  obstacles  to  create  a  violent  system  of  shock  waves  and  flames.  This  complex  is  a  precursor 
to  the  formation  of  a  detonation  through  a  process  called  DDT,  the  deflagration-to-detonation  transition. 
The  most  violent  type  of  a  steady  combustion  wave,  a  detonation,  consists  of  interacting  shock  waves  whose 
velocities  and  structure  are  determined  by  closely  coupled  and  intense  chemical  energy  release.  Detonations 
are  significantly  more  powerful  than  flame-generated  shocks  that  occur  earlier  before  the  transition.  In 
summary,  the  result  of  a  simple  spark  can  be  a  violent  deflagration  driving  shocks  or  even  more  violent 
detonation  waves.  Impacts  from  inert  shocks  produce  high  pressures  on  surrounding  walls,  and  even  higher 
pressures  can  be  generated  by  detonation  impact. 

Thus  a  simple  spark  can,  in  principle,  lead  to  a  series  of  combustion  waves  and  transitions  among  them 
that  result  in  enormous  stress  on  internal  mine  structures,  such  as  the  protective  seals  designed  to  contain 
natural  gas  in  unused  portions  of  the  mine.  For  this  reason,  the  possible  occurrence  of  DDT  and  a  detonation 
in  a  coal  mine  must  be  carefully  and  scientifically  evaluated  and  understood.  The  DDT  process  itself  must  be 
understood.  We  will  need  to  readjust  the  current  small-scale  notions  of  methane  detonability  to  large-scale 
actualities  of  confined  NG  explosions. 
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Most  recently,  NIOSH  IC9500  [1]  examined  the  fundamental  question:  Can  high  pressures,  DDT,  or 
detonations  develop  in  various  NG-air  mixtures,  and  what  factors  control  this  development?  This  report 
presented  results  of  thermodynamic  calculations  performed  to  establish  bounds  on  explosion  pressures  that 
can  develop  from  ignition  of  an  explosive  NG-air  mixture  in  mine  tunnels.  It  also  used  empirical  rules 
to  establish  conservative  estimates  of  the  run-up  length  or  distance  from  ignition  point  to  DDT  for  an 
explosion  in  a  mine  tunnel.  The  report  also  described  basic  gas-explosion  modeling,  using  commercially 
available  software,  which  provided  information  on  the  pressure-time  history  that  could  develop  from  a  NG- 
air  explosion.  These  models,  however,  did  not  contain  enough  detail  to  capture  the  dynamical  range  of  the 
DDT  process  or  predict  detonations,  nor  could  they  correctly  capture  the  liighest-pressure  blast  or  shock 
waves.  Nevertheless,  they  were  correct  in  indicating  that  very  high  pressures  have  developed,  and  this  was 
consistent  with  the  thermodynamic  calculations.  Another  group  of  studies  calculated  explosion  pressures 
created  by  increasingly  large  explosions  of  mixtures  of  NG  and  air  that  could  form  directly  behind  a  leaking 
seal.  None  of  the  codes  were  designed  to  be  able  to  predict  DDT.  They  provided  guidance  for  monitoring  and 
establishing  the  allowable  amount  of  explosive  mix  that  can  exist  behind  a  seal  of  given  strength,  without 
taking  into  account  the  possibility  of  detonations. 

NIOSH  IC9500  led  to  several  critical  observations  and  recommendations.  In  summary: 

1.  A  detailed  picture  of  the  temporal  and  spatial  evolution  of  the  pressure  within  a  typical  mining 
environment  did  not  exist,  and  the  situation  was  particularly  true  for  environments  in  which  a 
transition  from  a  deflagration  to  a  detonation  (DDT)  can  occur. 

2.  DDT  was  not  yet  fully  understood  for  conditions  typically  encountered  in  mine  tunnels. 

3.  It  is  important  to  take  action  to  understand  how  detonations  may  arise,  quantify  their  dangers 
as  accurately  as  possible,  and  use  this  information  to  determine  guidelines  for  protecting  mining 
personnel  and  the  mines  themselves. 

4.  We  would  like  to  be  able  to  detect  the  danger  through  monitoring  before  an  explosion  occurs  and 
take  action  to  mitigate  this  threat. 

Based  on  these  issues  and  observations,  the  Laboratory  for  Computational  Physics  and  Fluid  Dynamics 
(LCP&FD)  at  the  US  Naval  Research  Laboratory  (NRL),  in  close  coordination  with  the  NIOSH,  OMSHR 
(. National  Institute  for  Occupational  Safety  and  Health,  Office  of  Mine  Safety  and  Health  Research )  undertook 
an  effort  to  address  some  of  the  basic  questions  of  ignition,  flame  acceleration,  and  possible  DDT  in  scenarios 
relevant  to  coal  mines.  In  the  fifteen  years  prior  to  the  work  reported  here,  researchers  in  LCP&FD  developed 
a  series  of  numerical  simulation  tools  for  describing  the  behavior  of  gas-phase  explosions  and  the  transition 
to  detonation  from  first  principles.  Application  of  these  tools,  in  conjunction  with  carefully  performed 
laboratory  experiments,  has  both  validated  these  models  as  predictive  tools  and  created  a  body  of  new 
knowledge  on  the  behavior  of  explosions  in  gaseous  fuels.  An  important  component  for  success  in  this  earlier 
work  was  the  close  ties  maintained  between  experiments  and  simulations. 

Early  numerical  studies  using  the  LCP&FD  capabilities  showed  how  flames  in  low-pressure  mixtures  of 
acetylene  or  ethylene  and  air  can  accelerate  and  become  turbulent  due  to  repeated  shock  interactions,  and 
how  DDT  can  occur  in  fuel-air  systems.  Studies  were  also  performed  to  uncover  the  effects  of  boundary 
layers  and  wakes  on  flame  acceleration  and  possible  DDT.  All  of  these  early  simulations  were  based  on 
actual  laboratory  experiments  performed  at  the  detonation  research  facility  at  the  University  of  Wales  at 
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Aberystwyth,  UK,  and  they  had  the  remarkable  result  of  teaching  us  fundamental  mechanisms  and  explaining 
experimental  phenomena  that  were  not  previously  understood.  A  summary  review  of  the  work  done  up  until 
2003  is  given  in  [5]. 

Between  2003  and  2008,  the  LCP&FD  research  in  flame  acceleration  and  DDT  was  focused  on  simulating 
hydrogen  explosions,  specifically,  the  development  and  acceleration  of  flames  and  the  transition  to  detonations 
in  channels  with  obstructions.  The  purpose  of  this  work  was  to  produce  guidelines  for  safe  design  of  hydrogen 
storage  facilities.  A  partial  summary  of  the  results  of  this  work  is  given  in  [6-8].  An  important  aspect  of 
the  hydrogen-fuel  study  was  the  comparison  between  the  calculations  and  experiments,  which  were  carefully 
designed  to  describe  the  same  gases  and  have  the  same  geometrical  configurations.  The  numerical  simulations 
compared  well  with  experiments  and  provided  many  insights  into  flame  development  and  acceleration,  DDT, 
and  detonations  in  obstructed  channels.  The  geometries  considered,  an  example  of  which  is  given  in  Figure 
2.1,  are  relevant  to  the  problem  of  explosions  in  mine  tunnels  if  the  gas  were  changed  to  NG-air  and  the  size 
could  be  scaled  up  appropriately.  This  background  numerical  development  is  summarized  in  more  detail  in 
Section  2  of  this  report. 

In  the  three  years  of  the  current  project  with  NIOSH,  LCP&FD  personnel  have  focused  on  developing 
and  using  numerical  and  experimental  technologies  that  allow  us  to  address  critical  questions  related  to  NG 
explosions  in  coal  mines.  In  particular,  there  is  the  overriding  question: 

Given  a  large  enough  volume  of  a  flammable  mixture  of  NG  and  air ,  such  as  may  exist  in  a  coal 
mine,  can  a  weak  spark  ignition  develop  into  a  detonation? 

The  implications  of  an  answer  such  as  “yes,”  or  even  “sometimes,”  means  that  we  could  be  seriously  under¬ 
estimating  the  pressures  that  can  develop  after  an  explosion  in  a  mine  tunnel.  To  address  this  question  and 
its  implications,  LCP&FD  has  worked  on  two  major  tasks,  one  involving  numerical  simulation  and  the  other 
in  supporting  the  design  and  implementation  of  a  new  experimental  test  facility. 

The  work  in  numerical  simulation  has  focused  on  the  development  and  evolution  of  explosions  in  scenarios 
appropriate  to  those  in  mine  tunnels.  We  consider  methane,  a  primary  component  of  NG,  which  itself  can 
vary  in  composition  and  contain  anywhere  from  0.5-8%  of  higher  hydrocarbons  (HHC).  A  few  percent  of 
HHC  can  make  NG  easier  to  ignite  and  explode.  At  this  stage  of  the  model  development  and  testing,  however, 
we  are  considering  NG  and  methane  to  be  the  same.  In  this  work,  we  often  refer  to  NG-air  mixtures  as 
methane-air  mixtures.  Then  the  specific  tasks  included  in  this  work  are: 

1.  Develop  and  test  a  model  for  stoichiometric  methane-air  combustion.  Perform  a  series  of  simulations 
of  flame  propagation,  acceleration,  and  possible  DDT  in  channels  with  obstacles. 

2.  Develop  and  test  similar  models  for  fuel-rich  and  fuel-lean  combustion. 

3.  Examine  the  effects  of  spatially  varying  fuel  concentrations  in  a  simulation. 

4.  Document  the  effects  of  stochasticity,  i.e.,  examine  how  large-scale  uncertainty  in  system  behavior 
caused  by  minute  perturbations  in  initial  conditions  or  numerical  procedures,  affect  flame  acceleration 
and  DDT. 

Thus,  as  part  of  this  effort,  we  developed  appropriately  simplified  reduced  chemical-diffusion  models  that 
are  necessary  for  large-scale  multidimensional  unsteady  numerical  simulations.  These  models  were  then  used 
to  perform  benchmark  computations  of  the  flame  evolution  to  supersonic  deflagrations,  and,  if  appropriate, 
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detonations  in  methane-air  mixtures.  Computational  results  were  compared  with  pre-existing  experimental 
data.  Results  and  interpretations  of  these  simulations  are  described  in  detail  in  a  number  of  published  works 
[9-13]  and  in  subsequent  sections  of  this  report. 

The  experimental  effort  required  close  coordination  between  LCP&FD  personnel  and  NIOSH  engineers 
in  order  to  design  and  use  the  new  explosion  tube  at  Lake  Lynn  Laboratories  (LLL).  This  is  currently  the 
largest  operational  test  facility  for  NG  explosions  [14].  Its  successful  construction  required  design  of  new 
mixing,  diagnostic,  and  ignition  systems.  Its  use,  to  date,  has  lead  to  important,  new  results  concerning  the 
detonability  of  large  volumes  of  NG-air  mixtures.  More  details  of  the  operation,  tests,  and  implications  of 
these  results  will  be  discussed  in  Section  4  of  this  report. 

We  have  shown  that  the  complexity  and  scale  of  a  NG  explosion  preclude  us  from  having  any  analytic  or 
simplified  theoretical  analysis  of  the  complete  explosion  process,  from  initiation  through  to  possible  DDT.  It  is 
also  difficult,  dangerous,  and  extremely  expensive  to  study  large-scale  explosions  in  a  controlled,  experimental 
environment.  Because  of  the  practical  consequences  of  both  unwanted  and  purposeful  explosions,  we  must 
take  advantage  of  every  available  tool  we  have  for  studying  and  explaining  explosions.  We  must  combine 
numerical  simulations  with  carefully  designed  experiments  and  analyses,  both  to  understand  and  then  to 
use  what  we  have  learned  about  explosions  to  design  safe  and  efficient  systems.  This  report  summarizes  the 
work  done  by  NRL  and  OMSHR  to  date  to  create  a  reliable,  predictive  capability. 
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Figure  1.1.  Computed  flame  evolution  in  an  obstructed  channel  filled  with  a  stoichiometric  methane-air 
mixture  [9] .  Each  frame  is  a  window  102  long  cm  by  26  cm  high.  The  window  moves  through  the  system  to 
follow  the  reaction  front  and  the  flames  and  shocks  directly  behind  it.  The  height  of  the  channel  is  d  =  52  cm, 
so  that  the  half-channel  height  used  to  define  the  computational  domain  of  the  simulation  is  d/2  =  26  cm. 
The  blockage  ratio  is  defined  as  BR  =  2 h/d,  where  h  is  the  height  of  an  obstacle  and  d  is  the  heigh  of  the 
channel.  Here  BR  =  0.3. 
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Figure  1.2.  This  is  an  illustration  taken  from  a  computation  of  the  deflagration-to-detonation  transition, 
or  DDT,  occurring  in  array  of  obstacles  surrounded  by  a  stoichiometric  mixture  of  hydrogen  and  oxygen  at 
atmospheric  temperature  and  pressure  [8].  Each  obstacle  is  1  cm  in  each  direction.  The  background  is  in 
dark  blue.  The  leading  shocks  (light  blue)  are  followed  by  a  turbulent  flame.  See  Figure  2.6  for  more  details. 
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2.  Background:  Computing  Explosions 


Two  of  the  fundamental  questions  of  numerical  simulation  are:  How  do  we  solve  the  governing  equations? 
How  well  are  we  solving  them?  We  generally  agree  that  descriptions  of  explosions  in  normal  atmospheres 
require  a  solution  of  the  unsteady  Navier-Stokes  equations  (NSE)  with  some  representation  of  the  energy- 
release  process.  We  refer  to  these  coupled  fluid-chemical  equations  as  the  reactive  Navier-Stokes  equations 
(RNSE),  which  we  write  here  as 

g+V-(pv)=  0,  (2.1) 


d(pv) 

dt 


+  V  •  (pvv)  +  VP  +  V  •  r  =  0  , 


—  +  V  •  ({E  +  P) U)  +  V  •  (v  •  r)  +  V  •  (K  VP)  +  p q  Cl  =  0  , 


(2.2) 

(2.3) 


d(pY) 

dt 


+  V  •  (pYv) 


V-(pPVF)-pH  =  0  , 


(2.4) 


These  equations  describe  the  conservation  of  total  mass,  momentum,  and  energy,  and  the  balance  of  chemical 
species.  Here  p  is  the  mass  density,  v  is  the  fluid  velocity,  P  is  the  pressure,  E  is  the  energy  density,  Cl  is 
the  reaction  rate,  q  is  the  total  chemical  energy  release,  K  is  the  thermal  conduction  coefficient,  Y  is  mass 
fraction  of  fuel,  D  is  the  mass  diffusion  coefficient,  and  t  is  time.  The  viscous  stress  tensor  r  is  defined  as 


r  = 


pu  (V  •  v)  I  -  (Vv)  -  (Vv)tj 


(2.5) 


where  v  is  the  kinematic  shear  viscosity,  I  is  a  unit  matrix,  and  superscript  f  indicates  matrix  transposition. 

To  close  this  set  of  equations,  we  must  add  the  equation  of  state  and  the  definition  of  Cl.  The  chemical 
reaction  rate  can  be  described  by  a  complex  set  of  coupled  partial  differential  equations  representing  the 
conversion  of  reactants  into  products  and  the  energy  release,  or  it  can  be  approximated  by  a  single-step 
model  with  only  one  equation. 

What  differs  among  the  various  explosion  codes  is  complexity  of  the  chemical  model  and  the  equation 
of  state,  the  method  of  solution  of  RNSE,  how  this  method  is  implemented  in  the  code  and  computer,  how 
different  terms  in  the  equations  are  coupled,  and  how  complicated  geometry  and  numerical  resolution  issues 
are  handled.  The  requirements  for  solution  of  RNSE  in  various  regimes  of  explosions  are  described  below  in 
more  detail. 

The  material  presented  in  this  section  is  divided  into  three  parts.  First,  we  give  the  necessary  background 
for  understanding  the  evolution  and  development  of  numerical  methods  that  are  applicable  for  simulating 
explosions.  This  will  also  introduce  some  of  the  necessary  computational  terminology  and  concepts  used 
subsequently  in  this  report.  Second,  we  describe  three  numerical  models  used  previously  to  simulate  coal 
mine  explosions,  and  discuss  their  capabilities  and  limitations  for  computing  the  full  range  of  explosion 
stages.  Finally,  we  summarize  the  NRL  background  in  developing  and  using  models  for  numerical  simulation 
of  high-speed  reactive  flows. 
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2.1.  Overview:  Evolution  of  Numerical  Methods 


The  importance  of  having  the  ability  to  use  numerical  methods  to  solve  an  explosion  problem  was  realized 
in  the  mid  1900’s  when  understanding  the  complexities  of  an  actual  explosion  became  an  important  issue 
in  weapons  design  and  analysis.  These  early  efforts  produced  the  first  fmite-diference  methods  for  solving 
the  NSE.  All  of  the  these  early  methods  had  intrinsic  difficulties  generally  with  achieving  sufficient  accuracy, 
and  with  numerical  stability  specifically  when  they  were  used  to  compute  the  properties  of  shock  waves  and 
detonations  [15,  16]. 

The  intrinsic  difficulties  were  based  on  the  fact  that  the  finite-differencing  methods  had  an  inherent  nu¬ 
merical  instability  that  introduced  essentially  uncontrolled  numerical  oscillations  at  the  shock  discontinuity. 
The  first  attempt  to  fix  this  by  adding  artificial  damping  terms  (numerical  diffusion)  did  not  always  work 
as  expected  and  required  ad  hoc  parameters  to  control  the  added  numerical  diffusion  at  the  shock  front  [17, 
18].  This  method  also  caused  more  numerical  problems  when  it  was  coupled  to  highly  exothermic  chemical 
reactions  that  might  occur  near  shock  fronts  in  detonation  waves. 

The  major  algorithmic  breakthrough  in  our  ability  to  solve  the  fluid-dynamic  issues  with  computing 
explosions  occurred  in  the  early  1970’s  with  the  invention  of  flux- limiting  methods  [19].  These  methods 
addressed  the  numerical  origin  of  the  uncontrolled  numerical  diffusion  around  propagating  discontinuities 
by  imposing  solution  constraints  based  on  maintaining  conservation,  monotonicity,  and  causality.  Over  the 
next  ten  years,  many  communities  adopted,  expanded,  and  improved  the  new  approach.  Now  all  modern 
shock-capturing  methods  for  solving  the  NSE  with  shocks  either  directly  or  indirectly  use  flux  correction, 
especially  near  fluid-dynamic  discontinuities.  Application  of  these  methods  to  astropliysical  problems  in 
addition  to  weapons  research  further  aided  their  development  [20]. 

With  advances  in  computer  capabilities  (such  as  speed,  memory,  and  architecture),  the  ability  to  perform 
meaningful  numerical  simulations  advanced  enormously.  When  the  new  methods  for  solving  the  NSE  were 
combined  with  new  methods  for  representing  chemical  reactions  and  energy  release,  applications  opened  up 
in  almost  every  area  of  science,  including  atmospheric  and  solar  physics,  combustion,  and  astrophysics,  and 
almost  every  branch  of  engineering.  The  accuracy  of  numerical  solutions  became  mostly  limited  by  speed 
and  memory  of  available  computers. 

The  next  major  step  in  computation  was  the  recognition  that  when  a  fluid  dynamics  problem  is  solved 
with  a  monotone  method  (i.e.,  one  that  is  designed  so  that  it  does  not  create  any  new,  unphysical  maxima  or 
minima  in  the  solution) ,  there  is  often  no  need  to  resolve  the  viscous  scale  or  add  special  models  to  describe 
the  effects  of  turbulence.  This  concept  is  called  MILES  (monotone-integrated  large-eddy  simulation)  or  ILES 
(implicit  large-eddy  simulation)  [21],  This  is  based  on  an  interesting  property  of  monotone  methods:  at  the 
computational  grid  scale,  monotone  methods  dissipate  energy  in  a  very  physical  way  that  mimics  the  natural 
viscous  dissipation.  To  model  a  turbulent  flow,  it  is  sufficient  to  resolve  some  portion  of  the  inertial  range. 
The  flow  is  then  resolved  from  the  largest  scale  down  to  the  grid  scale,  at  which  the  energy  is  dissipated  by 
a  numerically  imposed  viscosity.  Thus  the  viscous  dissipation  scale  is  effectively  replaced  by  a  grid  scale. 

The  use  of  MILES  methods  means  that  the  extreme  turbulence  in  an  explosion  can,  in  principle,  be 
accurately  computed  without  resorting  to  ad  hoc  models  with  empirical  constants.  This  concept  and  its 
proof  have  taken  many  years  of  testing  to  validate,  but  now  we  believe  we  have  a  good  understanding  of  how 
and  why  it  works  [22] .  It  is  now  widely  accepted  that  if  we  resolve  enough  scales  in  a  turbulent  flow,  we  do 


not  have  to  resolve  all  length  scales  down  to  the  viscous  dissipation  scale. 

The  final  numerical  advance  that  has  changed  our  ability  to  compute  explosions  is  the  advent  of  adaptive 
gridding  [23,  24].  Explosion  problems  are  often  characterized  by  a  large  disparity  in  physical  scales  and  a 
concentration  of  important  small-scale  features,  such  as  flames,  shocks,  and  detonation  waves,  in  relatively 
small  areas  of  the  computational  domain.  Since  small-scale  features  require  small  computational  cells,  a 
uniformly  spaced  computational  mesh  throughout  the  entire  computational  domain  usually  means  that  there 
are  large  volumes  of  space  with  many  small  grid  cells  and  nothing  much  happening.  To  increase  the  efficiency 
of  computations,  it  makes  sense  to  use  fine  grid  around  important  flow  features,  and  coarse  grid  everywhere 
else.  Because  of  the  dynamic  nature  of  explosions,  this  adaptive  gridding  needs  to  be  dynamically  adjusted 
to  resolve  evolving  flow  features.  This  became  possible  with  developments  of  special  algorithms,  generally 
called  methods  for  adaptive  mesh  refinement  (AMR). 

As  we  show  below,  using  explosion  codes  based  on  monotone  methods  with  adaptive  mesh  refinement, 
coupled  with  the  advances  in  computer  power,  has  allowed  enormous  improvement  in  our  ability  to  compute 
an  explosion. 

2.2.  Previous  Simulation  of  Coal-Mine  Explosions 

There  are  three  main  requirements  for  an  RNSE  code  that  is  used  to  describe  the  full  range  of  states  that 
can  occur  in  NG  explosions. 

1.  The  numerical  model  must  solve  the  RNSE  using  a  method  appropriate  for  and  accurate  enough  to 
capture  the  development  and  evolution  of  turbulence  and  shock  phenomena. 

2.  There  must  be  an  “accurate  enough”  model  for  energy  release,  which,  in  our  case,  would  represent 
natural-gas  combustion. 

3.  The  computational  domain  and  what  it  contains  must  be  specified  and  resolved  numerically.  This  means 
that  the  code  needs  to  be  given  boundary  and  initial  conditions  that  make  the  problem  unique  and  apply 
to  a  specific  scenario. 

Many  adequate  numerical  methods  exist  for  solving  the  RNSE,  and  many  of  these  have  been  applied  to 
NG  explosions.  If  they  are  to  predict  DDT,  however,  they  need  to  be  monotone,  have  a  chemical  model  that 
can  reproduce  properties  of  both  laminar  flames  and  detonations,  and  have  enough  resolution  to  compute 
the  transition  from  laminar  to  turbulent  flames  to  detonations,  as  it  occurs.  One  of  the  main  restrictions  of 
many  of  the  current  commercial  explosion  codes  is  that  they  use  “turbulence  models.”  The  underpinning 
idea  of  a  turbulence  model  is  to  solve  the  mean,  large-scale  fluid-dynamics  equations  on  a  computational 
grid  in  which  the  smallest  grid  size  is  substantially  larger  than  the  scale  on  which  viscous  dissipation  occurs. 
Then  the  effects  of  sub-grid-scale  turbulence  and  viscous  dissipation  are  modeled  by  computing  extra  terms 
defined  by  the  turbulence  model  (see,  for  example,  [25]).  The  intended  effect  is  that  the  turbulence  model 
should  modify  the  RNSE  to  account  for  unresolved  small-scale  effects,  but  in  many  cases  this  adds  significant 
amounts  of  numerical  diffusion  that  can  mask  DDT.  Another  issue  related  to  the  use  of  turbulence  models, 
is  that  capturing  DDT  requires  resolving  the  problem  over  a  wide  range  of  time  and  space  scales.  This  could 
be  as  extreme  as  spatial  scales  ranging  from  micrometers  or  millimeters  (e.g.,  viscous  dissipation  scales, 
flame  thicknesses)  to  kilometers  (channel  lengths)  and  time  scales  ranging  from  nanoseconds  to  seconds 
(e.g.,  chemical  reaction  times). 
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An  excellent  discussion  of  many  of  the  gas  explosion  models  and  codes,  with  an  evaluation  of  their 
capabilities  and  limitations,  is  given  in  reviews  by  Lea  and  Ledin  [26]  and  Popat  et  al.  [27].  Commercial 
codes,  such  as  FLACS  [28,  29]  and  AutoReaGas  [30],  have  been  used  for  at  least  fifteen  years  by  a  number  of 
communities.  SAGE  [31]  is  a  code  containing  some  higher-level  algorithms  and  flexible  gridding  capabilities. 
All  of  these  codes  have  been  used  to  model  coal  mine  explosions  and  help  to  explain  observed  damage 
[1,  4,  32].  They  use  different  underlying  numerical  methods  and  implement  somewhat  different  chemical 
submodels.  Here  we  briefly  describe  relevant  features  of  these  codes  so  that  they  can  be  compared  and 
contrasted  to  each  other  and  to  the  the  NRL  codes  and  models  we  are  currently  using.  This  comparison  will 
help  to  explain  why  we  need  a  different  code  and  different  approaches  to  capture  DDT. 

Since  both  AutoReaGas  and  FLACS  are  based  on  monotone  fluid-dynamics  solvers,  they  should,  in 
principle,  be  able  to  compute  the  properties  of  explosions  and  shocks.  AutoReaGas  uses  a  high-order 
monotone  algorithm  (based  on  a  version  of  the  flux-corrected  transport  (FCT)  algorithm  [16,  19]  developed 
at  NRL)  and  FLACS  uses  a  lower-order  algorithm.  Practically,  the  higher-order  algorithm  in  AutoReaGas 
means  that  this  code  requires  fewer  computational  cells  to  resolve  a  feature  in  the  flow,  but  it  also  means 
that  it  could  take  slightly  longer  to  perform  the  calculation.  The  turbulence  models  used  in  both  codes  are 
the  k  —  e  models  [25],  which  relate  the  dissipation  rate  of  unresolved  turbulence  kinetic  energy  ( k )  to  the 
production  of  turbulence  kinetic  energy  in  the  modified  form  of  the  RNSE.  Based  on  the  previous  discussions, 
a  question  arises:  Why  does  a  monotone  method  need  a  turbulence  model?  The  answer  must  be  that  the 
computations  performed  with  the  code  are  not  resolved  enough  to  take  advantage  of  the  natural  dissipation 
properties  of  the  algorithm  they  are  using.  This  is  not  helped  by  the  use  of  additional  diffusion  terms  that 
come  from  the  turbulence  models. 

In  the  AutoReaGas  or  FLACS  calculations,  the  smallest  grid  size  is  of  the  order  of  0.5  m  and  is  much 
larger  than  a  flame  thickness,  which  is  about  1  mm  for  a  stoichiometric  methane-air  mixture  [33].  Thus  part 
of  the  turbulence  model  must  represent  features  of  chemical  reactions  and  energy  release  that  take  place  on 
the  scale  of  the  laminar  flame  thickness.  The  resulting  description  of  the  turbulent  flame  itself  will  be  some 
averaged  profile,  but  there  are  constants  in  the  turbulence  and  chemical  models  that  help  ensure  that  the 
correct  amount  of  energy  is  released  overall,  and  this  ensures  the  correct  flame  temperatures  and  deflagration 
velocities. 

AutoReaGas  and  FLACS  are  RNSE  codes  because  they  contain  combustion  submodels  that  describe 
the  concentration  change  rates  of  reactant  and  product  species  and  the  associated  energy  release  rate.  For 
large-scale  explosion  simulations,  these  submodels  are  empirical.  For  instance,  AutoReaGas  uses  an  empirical 
correlation  between  reaction  rate  and  flame  speed.  FLACS  uses  the  “/?  flame  model”  that  correlates  turbulent 
burning  velocity  with  turbulence  parameters.  In  both  models,  an  increase  in  turbulent  kinetic  energy  results 
in  an  increase  in  the  reaction  rate.  These  models,  however,  do  not  properly  consider  the  physics  of  detonation 
or  DDT. 

AutoReaGas  and  FLACS  are  designed  to  approximately  model  deflagrations  and  can  be  applied  to 
relatively  low-pressure  explosions.  In  most  applications  of  these  codes  in  the  oil,  gas,  and  chemical  industries, 
the  computed  and  measured  explosion  pressures  do  not  exceed  about  0.5  MPa.  The  large  computational 
cell  sizes  they  use  mean  that  details  of  turbulence  and  its  effect  on  flame  acceleration  and  explosion  pressure 
are  not  computed.  The  coarse  grid  also  tends  to  “smear”  or  diffuse  shocks  and  other  rapidly  rising  pressure 
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waves  that  are  important  for  DDT  mechanisms.  Thus,  these  models  cannot  predict  extremely  high  pressures 
that  could  occur  in  a  mining  explosion  as  a  result  of  possible  DDT. 

Despite  these  shortcomings  at  high  pressures,  such  models  still  provide  useful  insights  into  many  prac¬ 
tical  applications  of  interest  at  lower  pressures.  FLACS  and  AutoReaGas  have  been  used  to  help  analyze 
explosions  in  industrial  facilities,  such  as  oil  production  platforms,  oil  refineries,  chemical  process  plants,  and 
mines.  In  particular,  these  codes  have  been  used  to  simulate  explosions  in  off-shore  oil  and  gas  production 
facilities  since  the  Piper- Alpha  disaster  in  1988  [34,  35].  Several  groups  in  Europe  have  made  attempts  to 
use  these  models  to  study  gas  explosions  in  mines.  NIOSH  has  used  them  with  success  to  help  analyze  mine 
explosions  in  LLEM,  which  probably  represents  the  most  extensive  use  of  these  models  in  a  mining  industry 
application.  The  full  extent  of  this  initial  gas  explosion  modeling  application  in  mining  is  well  documented 
in  several  reports  [36,  37]. 

SAGE  code  [31]  is  used  by  DOD  to  study  blast  effects  from  weapons.  It  uses  a  relatively  high-order 
monotone  method  for  capturing  shocks,  and  to  this  it  adds  a  method  for  adaptive  mesh  refinement  (AMR) 
so  that  it  can  capture  developing  features  in  the  flow.  It  also  contains  models  that  allow  it  to  evaluate 
the  effects  of  explosions  on  solid-fluid  interfaces,  such  as  between  concrete  and  air.  It  can  simulate  strong 
shock  waves  and  complex  combustion  fronts  because  AMR  can  introduce  resolution  where  it  is  needed.  It 
does  not  use  a  turbulence  model,  but  again,  because  it  is  based  on  a  monotone  algorithm  with  good  grid 
resolution,  it  can  take  advantage  of  MILES.  The  chemistry  model  is  based  on  an  Arrhenius  reaction  model 
that  describes  the  combustion  rate  in  a  computational  cell  in  terms  of  the  concentrations  of  reactants  and 
products  and  chemical  rate  constants.  The  Griineisen  equation-of-state  relates  pressure  to  the  density  and 
the  internal  energy  of  the  combustion  products  for  methane-air  burning.  Because  of  its  versatility,  SAGE  can 
model  combustion  through  either  a  deflagration  or  detonation  mechanism.  It  is  not  known  if  it  can  compute 
DDT  from  first  principles,  or  whether  the  chemical  model  is  adequate  to  do  this.  It  is  also  not  known  if  the 
type  of  adaptive  gridding  it  uses  would  be  adequate  for  computing  DDT.  Currently,  the  transition  from  a 
deflagration  to  a  detonation  is  done  in  an  ad  hoc  manner. 

It  is  probably  accurate  to  say  that  complete  reactive  CFD  simulation  of  any  realistic  large-scale  explo¬ 
sion  is  beyond  any  of  today’s  computer  capabilities,  and  probably  beyond  our  ability  to  specify  the  initial 
conditions  and  input  constants.  The  only  viable  approach  is  to  select  levels  of  algorithms,  models,  and 
gridding  so  that  those  chemical  and  physical  processes  that  we  think  are  important  are  resolved  and  the 
combustion  phenomena  they  create  are  captured. 

2.3.  Previous  NRL  Explosion  Simulations:  Deflagrations,  Detonations,  and  DDT 

The  computational  tasks  performed  for  NIOSH  rest  heavily  on  prior  work  at  NRL  in  algorithm  development, 
code  construction,  and  use  of  these  codes  on  high-performance  computers  to  study  fundamental  mechanisms 
of  explosions.  The  earliest  NRL  work  includes  some  of  the  first  computational  studies  of  propagating  flames 
and  detonations,  and  culminates  in  the  publication  of  the  book  Numerical  Simulation  of  Reactive  Flow  [15]. 
Since  ^1995,  there  has  been  a  dedicated  effort  to  understand  DDT.  The  parts  of  that  work  which  are  relevant 
to  this  report  comprise  the  rest  of  this  Section. 

One  general  conclusion,  however,  is  that  it  is  possible  to  compute  the  DDT  process  by  solving  the  RNSE, 
given  enough  numerical  resolution  of  critical  regions  in  the  system  and  a  chemical  model  that  spans  the 
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combustion  regimes  from  autoignition  to  laminar-flame  propagation.  We  come  back  to  this  major  conclusion 
at  the  end  of  this  section. 

2.3.1 .  Early  Development  of  Codes:  Simulations  of  Flames  and  Detonations 

The  earliest  work  at  NRL  in  numerical  simulations  of  reactive  flows,  and  the  particular  application  to  ex¬ 
plosions,  began  at  the  end  of  the  1960’s  with  simulations  of  high-altitude  nuclear  bursts.  This  required 
constructing  numerical,  physical,  and  nonequilibrium  chemical  models  for  effects  of  blasts  in  air.  The  ion¬ 
ization,  radiation,  and  products  of  the  chemical  reactions  produced  the  signatures  of  the  explosion.  For  our 
current  concerns  and  discussions  of  DDT,  the  significant  advances  that  came  from  this  early  work  are  the 
development  of  flux-limiting  algorithms  and  their  application  for  the  NSE,  fast  solution  methods  for  ordinary 
differential  equations  to  represent  the  chemical  reactions,  and  the  new  numerical  algorithms  that  allowed 
simulation  of  high-speed  reactive  flows. 

In  the  mid-1970’s,  NRL  work  turned  to  combustion,  which  emphasized  chemically  reactive  flows  where 
hydrogen  or  hydrocarbons  were  the  fuels  consumed.  Much  of  this  work  is  summarized  in  [15],  which  describes 
methods  and  applications.  Subsequently,  many  of  these  approaches  were  developed  further  at  NRL  and 
elsewhere  to  the  point  where  they  had  been  incorporated  into  commercial  CFD  codes  and  used  in  industrial 
situations.  The  second  edition  of  this  book,  which  was  published  in  2001  [16],  also  gave  information  on  how 
to  use  these  commercial  and  other  acquired  codes.  NRL  work  in  explosions,  deflagrations,  and  detonations 
focused  on  computing  the  single-  and  nrulti-dimensional  structure  of  propagating  flames  and  detonations, 
with  some  consideration  of  initiation  and  quenching.  In  addition,  there  was  related  work  on  chemical-fluid 
instabilities,  such  as  those  which  led  to  chemical-acoustic  coupling  [38] . 

The  following  summary  of  work  on  DDT,  starting  around  1995  and  continuing  through  to  the  present, 
is  cast  in  terms  of  major  conclusions  with  supporting  results  shown  in  figures  with  detailed  captions.  Only 
those  aspects  of  the  conclusions  and  observations  that  are  related  directly  to  DDT  in  methane-air  and  natural 
gas-air  mixtures  will  be  discussed  here. 

2.3.2.  Fundamentals  of  Deflagration-to-Detonation  Transition 

The  overriding  objective  of  the  simulations  was  to  determine  if  we  could  compute  DDT  from  “first  principles.” 
More  specifically,  we  wanted  to  determine  if  we  could  resolve  a  large  enough  range  of  the  temporal  and  spatial 
scales  to  see  the  transition.  At  the  minimum,  this  would  require  resolving  spatial  scales  from  the  laminar 
flame  thickness  to  the  size  of  the  system,  when  and  where  they  are  important  in  the  process.  The  approach 
we  used  combined  numerical  simulations  with  comparisons  to  laboratory  experiments  in  an  attempt  to 
find  mechanisms  for  unexplained  observations.  The  primary  tool  was  multidimensional,  time-dependent 
numerical  simulation  of  unsteady  combustion  using  high-performance  parallel  computers.  A  more  extensive 
description  is  given  in  [5],  and  presented  in  more  detail  in  the  references  therein. 

Conclusion:  We  can  compute  DDT  in  systems  with  shocks  and  flames  present. 

Figure  2.1  describes  one  of  the  first  problems  studied,  shock-flame  interactions  and  the  evolution  to  a 
turbulent  flame  and  DDT  in  an  enclosed  channel  filled  with  an  acetylene-air  mixture.  The  configuration  was 
based  on  experiments  performed  at  the  University  of  Wales  at  Aberystwyth.  This  result  proved  that  it  was 
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possible  to  simulate  DDT  by  solving  the  RNSE.  The  result  was  confirmed  by  convergence  studies  for  this 
flow  and  by  comparisons  to  experiments  [39].  The  mechanism  for  DDT  was  the  occurrence  of  a  hot  spot 
(or  ignition  center)  that  arose  in  shocked,  compressed,  and  heated  gas  in  the  vicinity  of  a  turbulent  flame. 
Subsequent  computations,  with  stronger  intensity  shocks  and  turbulence,  produced  hot  spots  closer  to  the 
flame  brush.  In  one  of  the  most  intense  cases,  DDT  occurred  in  funnels  of  unburned  material  inside  the 
flame  brush. 

Conclusion:  Detonations  arise  when  local  conditions  in  unreacted  material  allow  a  spontaneous  wave  to 
form,  and  this  wave  evolves  into  a  shock  that  is  strong  enough  to  become  a  detonation  which  can  propagate 
outside  the  gradient. 

The  gradient  mechanism  of  detonation  ignition,  proposed  by  Zeldovich  [40]  and  discussed  and  generalized 
by  many  researchers  since,  teaches  us  that  a  gradient  in  chemical  reactivity  can  generate  a  supersonic 
“spontaneous  wave”  that  may  undergo  transition  to  a  detonation.  If  it  does  not,  it  results  in  a  decoupled 
flame  and  a  shock.  Previous  theoretical  and  some  experimental  studies  of  this  ignition  mechanism  had  been 
done  in  isolation  from  a  turbulent  flame:  a  gradient  was  proposed  or  created  and  the  subsequent  behavior 
studied. 

Our  simulations  [41-43]  have  shown  that  such  gradients  arise  naturally  around  hot  spots  in  unreacted 
material  as  part  of  the  background  flow  in  which  the  turbulent  flame  exists.  Close  examination  of  the  hot 
spot  ignition  have  shown  that  it  occurred  in  a  gradient  of  reactivity  which  allowed  the  development  of  a 
spontaneous  wave. 

Conclusion:  Interactions  of  shocks  and  flames  are  important  in  creating  the  conditions  in  which  hot  spots 
can  arise  and  subsequently  produce  detonations. 

Flames  generate  shocks,  especially  in  confined  spaces  or  in  the  presence  of  obstacles.  In  turn,  shock 
interactions  with  flames  create  and  drive  the  turbulence  in  flames  far  more  dynamically  than  other  combustion 
instabilities.  Enhanced  turbulence  means  a  higher  energy-release  rate,  and  more  fluctuations  in  background 
material  that  promote  hot-spot  formation  and  increase  the  likelihood  of  DDT.  This  is  shown  in  Figure  2.1, 
and  discussed  in  more  detail  in  [41-43]. 

Conclusion:  Reflected  shocks  and  boundary  layers  can  interact  with  shocks  and  Barnes  to  reduce  the  time 
to  DDT. 

High-intensity  turbulent  deflagrations  in  channels  often  generate  shocks  and  high  speed  flows  that  in¬ 
teract  with  channel  walls  and  flames.  These  interactions  enhance  turbulence,  and  can  also  produce  complex 
dynamic  structures  that  allow  the  flame  to  spread  supersonically  relatively  to  the  unburned  material.  This 
detonationless  supersonic  flame  spread  significantly  increases  the  energy  release  in  the  system,  and  eventually 
accelerates  DDT. 

For  example,  a  shock  reflected  into  an  incoming  flow  can  interact  with  a  boundary  layer  and  produce  a 
complex  A-shock  structure  containing  a  leading  oblique  shock  followed  by  a  recirculation  region.  If  the  flame 
is  close  enough  to  this  structure,  it  becomes  entrained  in  the  recirculation  region  and  attached  to  the  oblique 
shock.  The  resulting  reactive  structure  appears  as  a  shock-flame  complex  propagating  at  approximately  one 
half  of  the  CJ  velocity.  This  leads  to  an  increase  in  the  energy-release  rate,  the  formation  of  Mach  stems  in 
the  middle  of  the  channel,  and  creation  of  multiple  hot  spots  behind  the  Mach  stem,  thus  facilitating  DDT 
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[44].  An  example  of  this  phenomenon  is  shown  in  Figure  2.2  that  compares  simulations  with  and  without 
the  effects  of  boundary  layers  for  a  shock-flame  interaction  in  an  etlrylene-air  mixture.  Similar  phenomena 
can  be  produced  by  interactions  of  reflected  shocks,  flames,  and  wakes  behind  obstacles  [45]. 

Conclusion:  For  the  large-scale  features  and  transitions  in  the  Bow,  2D  simulations  give  both  qualitative 
and  surprisingly  good  quantitative  agreement  with  3D  simulations  and  experiments. 

Figure  2.3  shows  the  results  of  a  calculation  with  boundary  layers  for  a  three-dimensional  configuration 
analogous  to  the  two-dimensional  version  shown  in  Figure  2.2.  For  this  ethylene-air  mixture,  2D  and  3D 
simulations  give  similar  times  to  DDT  [46].  A  large  part  of  this  agreement,  which  seemed  fortuitous,  was 
attributed  to  the  unusual  nature  of  the  nonequilibrium,  shock-generated  turbulence  in  the  flow  [44-46]. 

2.3.3.  Flame  Acceleration  and  DDT  in  Channels  with  Obstacles 

This  work  focused  on  flame  propagation,  acceleration,  and  eventual  DDT  in  obstructed  channels  filled  with 
a  stoichiometric  mixture  of  hydrogen  and  air.  The  geometrical  parameters  of  the  channel  and  the  method 
of  flame  ignition  were  varied  to  extract  general  principles  that  could  be  used  to  evaluate  the  likelihood  of 
DDT.  Issues  related  to  numerical  accuracy  were  addressed.  Validating  experiments  performed  at  Warsaw 
Institute  of  Technology  and  AberShock,  Inc.  (UK)  tested  selected  aspects  of  the  results  [47,  48].  Although 
the  hydrogen-air  mixture  used  in  this  work  is  much  more  sensitive  than  NG-air,  the  general  principles  learned 
can  be  used  to  guide  experiments  and  simulations  in  NG.  For  that  reason,  we  now  present  a  summary  of  the 
hydrogen  simulations. 

Conclusion:  A  spark  ignition  in  an  obstructed  channel  can  lead  to  DDT.  Detonations  can  form  when  shocks 
collide  with  obstacles. 

At  the  initial  stages  of  flame  propagation,  a  thermal  expansion  of  hot  combustion  products  generates  a 
flow  ahead  of  the  flame.  This  flow,  that  becomes  non-uniform  due  to  interactions  with  obstacles,  distorts  the 
flame,  thus  increasing  its  surface  area.  This  increases  the  energy  generation  rate  and  results  in  flame  and  flow 
acceleration.  Flame-vortex  interactions  and  a  number  of  fluid  instabilities,  such  as  Rayleigh-Taylor  (RT) 
and  Kelvin-Helmholtz  (KH),  continue  to  increase  the  flame  surface  area,  and  further  accelerate  the  flow. 
The  accelerating  flow  generates  shocks  that  reflect  from  obstacles  and  interact  with  flame  surface  causing 
Richtmyer-Meshkov  (RM)  instabilities  that  become  a  dominant  mechanism  for  the  flame  wrinkling.  The 
flame  continues  to  accelerate,  shocks  strengthen,  and  eventually  shock  reflections  from  obstacles  create  hot 
spots  that  produce  detonations  through  the  gradient  mechanism.  A  representative  example  of  this  process 
is  shown  in  Figure  2.4  [7]. 

A  series  of  simulations  performed  for  different  energies  of  ignition  sparks  shows  that  though  the  ignition 
energy  affects  the  initial  flame  development  and  changes  the  time  to  DDT,  this  often  has  little  or  no  effect 
on  the  distance  to  DDT  [7]. 

Conclusion:  As  the  height  of  the  channel  increases,  the  Bnal  combustion  regime  changes  from  choking 
deBagration  to  quasi-detonation,  and  then  to  detonation  [7]. 

For  narrow  channels,  the  choking  deflagration  is  the  final  regime  resulting  from  flame  acceleration.  In 
this  regime,  the  turbulent  flame  and  leading  shock  form  a  complex  in  which  they  remain  a  relatively  constant 
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distance  apart.  Together,  they  travel  at  about  (0.4-0. 5  Dcj).  The  leading  shock  diffracts  as  it  passes  over 
each  obstacle,  and  then  it  reflects  from  the  bottom  wall  after  each  diffraction. 

For  larger  channels,  the  choking  deflagration  with  the  shock-flame  complex  undergoes  a  transition  to 
a  quasi-detonation.  In  this  regime,  there  are  repeated  DDT  events  that  produce  detonations  as  shocks 
collide  with  obstacles.  These  detonations  do  not  survive  diffraction  over  an  obstacle,  but  they  are  constantly 
reignited  as  a  the  leading  shock  collides  with  the  next  obstacle. 

For  even  larger  channels,  the  DDT  leads  to  a  detonation  that  does  survive  diffraction  and  propagates 
through  the  system. 

Conclusion:  More  obstacles  per  unit  length  create  more  perturbations,  which  increase  the  flame  surface 
area  more  quickly.  Thus,  the  flame  speed  grows  faster. 

Conclusion:  Detonations  form  when  shocks  collide  with  obstacles.  These  collisions  are  most  effective  in 
generating  hot  spots  and  detonations  when  they  involve  Mach  stems. 

A  Mach  stem  may  form  when  the  diffracted  leading  shock  collides  with  the  channel  wall.  If  this  Mach 
stem  is  strong  enough,  its  collision  with  the  next  obstacle  may  produce  a  detonation.  The  strength  of  this 
Mach  stem  depends  on  the  strength  of  the  leading  shock,  the  height  of  the  obstacle  that  causes  the  shock 
diffraction,  and  the  distance  to  the  next  obstacle.  Decreasing  the  distance  between  obstacles  eventually 
eliminates  the  Mach  stem,  and  then  the  detonation  initiation  requires  a  stronger  leading  shock  [6]. 

Thus,  the  effects  of  obstacle  spacing  on  flame  acceleration  and  the  actual  detonation  initiation  are 
counteracting  to  some  extent.  This  is  why  increasing  the  number  of  obstacles  per  unit  length  may  increase 
the  distance  to  DDT,  even  though  the  flame  speed  would  grow  faster. 

Different  shock  configurations  that  form  just  before  the  detonation  initiation  are  illustrated  by  Figure 
2.5.  The  figure  shows  the  leading  shock  and  turbulent  flame  brush  developed  in  a  series  of  hydrogen-air 
simulations  in  which  the  spacing  between  successive  obstacles,  S,  inside  the  channel  is  systematically  varied. 
We  observe  two  distinct  configurations: 

Configuration  1:  For  S  equal  to  or  larger  than  the  channel  width  d  (S/d  =  1,  1.5,  2  in  Figure  2.5),  a  Mach 
stem  is  formed  by  the  diffracting  shock  reflecting  from  the  bottom  wall.  The  detonation  is  ignited  when  this 
Mach  stem  collides  with  an  obstacle. 

Configuration  2:  For  S  smaller  than  the  channel  width  (S/d  =0.25,  0.375,  0.5,  0.75  in  Figure  2.5),  Mach 
stems  do  not  form. 

The  detonation  initiation  for  this  configuration  is  more  difficult  since  the  leading  shock  without  a  Mach 
stem  may  not  be  strong  enough.  For  S/d  =  0.375,  0.5  and  0.75  detonations  do  not  develop.  Isolated 
detonations  sometimes  appear  in  pockets  of  unburned  material  between  obstacles  behind  the  leading  edge 
of  the  flame,  but  they  cannot  spread.  For  S/d  =  0.25,  however,  multiple  isolated  detonations  contribute  to 
the  final  stages  of  flame  and  shock  acceleration,  and  the  leading  shock  becomes  strong  enough  to  ignite  a 
detonation  by  a  direct  collision  with  the  top  part  of  an  obstacle.  The  distance  to  DDT  in  this  case  is  about 
3  times  longer  than  for  the  case  in  which  S/d  =  1. 
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2.3.4.  DDT  in  Large  Spaces  Containing  Obstacles 

Here  we  consider  a  large  area  with  many  obstacles,  either  square  or  rectangular,  with  different  kinds  of 
symmetrical  or  asymmetrical  placement.  This  is  work  in  progress,  so  that  any  conclusions  are  preliminary 
and  subject  to  revision  or  qualification  as  the  analysis  proceeds.  There  are,  however,  a  few  observations  that 
are  important  to  note  now  as  they  are  relevant  to  explosions  in  coal  mines. 

Conclusion:  The  evolution  from  a  spark  or  a  flame  to  a  detonation  in  a  two-dimensional  array  of  obstacles 
follows  an  overall  pattern  similar  to  what  was  seen  in  channels  with  obstacles. 

This  is  shown  by  the  example  in  Figure  2.6.  After  ignition,  a  small  flame  propagates  for  some  time 
in  a  laminar  regime,  and  then  accelerates  as  its  surface  bends  with  the  flow  around  obstacles.  Turbulence 
develops  as  shear  layers  appear  in  the  flow  and  become  unstable,  and  this  contributes  to  the  burning  rate 
increase.  As  the  flow  and  flame  accelerate,  shocks  are  formed  ahead  of  the  flame.  The  shocks  are  reflected 
by  obstacles  and  propagate  back  to  the  flame,  thus  causing  shock-flame  interactions  and  further  increasing 
the  energy-release  rate.  When  shocks  become  strong  enough,  their  collisions  with  obstacles  create  hot  spots 
that  produce  detonations  [8,  49]. 

Conclusion:  In  structured  arrays  of  obstacles,  there  are  directional  preferences  for  flame  spread  and  the 
occurrence  of  DDT,  and  these  directions  are  different  for  different  stages  of  flame  acceleration 

In  structured  arrays,  some  directions  are  more  obstructed  than  others,  and  the  flow  loses  more  momen¬ 
tum  in  more  obstructed  directions.  This  delays  the  flame  spread  in  more  obstructed  directions  at  initial 
stages  of  flame  acceleration.  On  the  other  hand,  the  number  of  shock-flame  interactions  is  higher  in  more 
obstructed  directions,  and  this  leads  to  a  higher  flame  acceleration  and  faster  strengthening  of  the  leading 
shocks  in  these  directions  at  later  stages.  As  a  result,  detonations  first  appear  in  more  obstructed  directions. 

Detonations  spread  in  all  directions  and  interact  with  obstacles  causing  shock  reflections  and  diffractions. 
Some  parts  of  detonation  fronts  are  quenched  by  diffractions,  but  they  are  reignited  by  collision  of  decoupled 
shocks  or  overtaken  by  other  detonations.  Thus  detonations  continue  to  spread  and  quickly  burn  all  the 
material  between  obstacles  [8]. 

2.3.5.  Computing  DDT:  General  Observation 

Conclusion:  Two-Dimensional  and  three-Dimensional  simulations  have  often  given  the  same  values  of  run¬ 
up  distances  to  DDT  consistent  with  experimental  results.  It  is  often  sufficient  to  perform  two-dimensional 
simulations. 

This  conclusion  is  based  on  repeated  comparisons  between  simulations  and  experiments,  and  between 
2D  and  3D  simulations  [5,  49].  It  appeared  early  in  our  studies,  and  has  been  repeatedly  borne  out  by 
new  comparisons.  The  physical  meaning  of  this  conclusion  is  a  complicated  concept.  Our  current  belief  is 
that  the  origin  of  this  agreement  lies  in  the  the  particular  fundamental  properties  of  the  nonequilibrium, 
non-Kolmogorov  turbulence  in  these  flows  that  are  somewhat  similar  in  2D  and  3D.  This  is  a  topic  which 
is  important  to  understand  and  is  currently  under  study  in  another  research  program.  Meanwhile,  we  are 
constantly  checking  the  assumptions  involved  and  use  this  property  to  our  advantage  in  the  simulations 
described  below  for  methane-air  mixtures, 
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Conclusion:  It  is  possible  to  compute  the  DDT  process  by  solving  the  RNSE,  given  enough  numerical 
resolution  of  critical  regions  in  the  system  and  a  chemical  model  that  spans  the  combustion  regimes  from 
autoignition  to  laminar-flame  propagation. 

Now  that  we  have  shown  this,  the  questions  are:  What  are  the  limitations?  How  can  we  beat  them? 
How  do  we  apply  the  information  we  have  and  could  obtain? 
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Figure  2.1.  This  figure  is  taken  from  one  of  the  first  computations  showing  DDT  arising  from  a  gradient 
of  reactivity  in  the  flow  [41].  The  schematic  on  the  top  shows  the  conditions  simulated,  which  were  modeled 
on  experiments  performed  at  the  University  of  Wales  Aberystwyth.  Time  sequence  of  density  fields  showing 
the  overall  evolution  of  a  Mach  1.5  incident  shock  interacting  with  a  low-pressure  acetylene-air  flame,  cor¬ 
responding  to  the  initial  and  boundary  conditions  shown  on  the  top.  The  computational  domain  for  each 
frame  is  32  cm  by  1  cm.  (The  vertical  scale  is  stretched  by  a  factor  of  1.28.)  Time  steps  are  shown  on  the 
left  side  and  physical  times  are  shown  on  the  right  side.  Note  that  time  differences  between  frames  are  not 
equally  spaced,  but  are  clustered  near  the  DDT  event  starting  at  step  2200.  Incident  shock,  I;  Reflected 
shocks  Rl,  R2,  R3. 

The  initial  frame,  step  0,  shows  the  shock  on  the  far  left  and  the  unperturbed  flame  approximately 
in  the  middle  of  the  domain.  As  time  passes,  this  incident  shock  I  approaches  the  flame,  and  the  shock- 
flame  interaction  begins  around  step  400.  When  I  interacts  with  the  left  side  of  the  flame,  the  result  is  a 
rarefaction  wave  that  moves  back  upstream,  and  a  transmitted  shock  that  moves  through  the  flame.  When 
this  transmitted  incident  shock  interacts  with  the  right  side  of  the  flame,  it  produces  a  weak  reflected  shock 
Rl  that  also  moves  upstream,  first  through  the  material  burned  by  the  flame,  and  then  through  the  intially 
shocked,  unburned  material.  The  rarefaction  wave  and  Rl  can  be  seen  in  frames  600  through  1700.  The 
incident  shock  transmitted  through  the  flame  moves  downstream  toward  the  end  wall,  as  seen  in  frames  600 
and  800,  and  reflects  from  the  wall  by  step  900.  This  produces  a  strong  reflected  shock  R2. 

The  interactions  of  I  and  Rl  with  the  flame  trigger  a  Richtmyer-Meshkov  (RM)  instability,  which  creates 
a  large  funnel  of  unburned  material  that  penetrates  the  burned  region.  At  frame  1000,  R2  begins  to  interact 
with  the  distorted  flame.  This  generates  a  new  sequence  of  reflected  and  transmitted  shocks  and  rarefactions 
which  further  distort  the  flame.  The  shock  R2  passes  through  the  flame  by  step  1200,  and  continues  to  move 
upstream.  It  finally  merges  with  Rl  between  steps  1700  and  1800.  When  R2  passes  through  the  left  side 
of  the  flame  at  step  1200,  it  generates  another  reflected  shock  that  moves  to  the  right.  This  shock  passes 
through  the  flame,  reaches  the  end  of  the  tube,  and  reflects  from  the  end  wall  between  steps  1300  and  1400, 
producing  the  R3.  The  shock  R3  passes  through  the  flame  again,  further  distorting  the  flame,  and  finally 
emerging  from  the  left  side  of  the  flame  at  step  1600.  This  third  emerging  shock  R3  propagates  in  material 
that  was  shocked  three  times  previously  (incident  shock  plus  two  reflected  shocks).  It  eventually  merges 
with  the  combined  reflected  shock  (R1-R2)  between  steps  1900  and  2000. 

Between  steps  2200  and  2260,  two  sudden  explosions  occur  almost  simultaneously  in  the  material  behind 
the  merged  reflected  shock,  R1-R2-R3.  These  explosions  coalesce  to  form  a  detonation  that  spreads  in  all 
directions.  As  it  moves  to  the  right,  it  enters  the  flame  between  steps  2400  and  2500,  and  decays  into  a 
shock. 
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Figure  2.2.  This  figure  is  taken  from  one  of  the  first  computations  of  the  effects  of  boundary  layers  on 
shock- flame  interactions  and  DDT  [44]  showing  DDT  arising  from  a  gradient  of  reactivity  in  the  flow  [41]. 
The  schematic  on  the  top  shows  the  conditions  simulated,  which  were  modeled  on  experiments  performed 
at  the  University  of  Wales  Aberystwyth.  The  sequence  of  density  fields  showing  shock-flame  interactions 
and  DDT  in  a  stoichiometric  ethylene-air  mixture  at  Po  =  1.44  x  103  Pa  for  initial  and  boundary  conditions 
shown  on  top.  Two  cases,  with  and  without  boundary  layer  were  computed.  Incident  shock  Mach  number 
Ms  =  1.9.  Time  (/zs)  is  given  on  the  left  side  of  each  frame.  The  letters  indicate  the  incident  shock  I,  the 
flame  F,  reflected  shocks  R1  and  R2,  funnels  of  unreacted  material  J1  and  J2,  detonations  D1  and  D2,  and 
bifurcated  structures  Bl,  B2  (not  shown)  and  B3. 

The  sequence  on  the  left  is  taken  from  a  simulation  with  free-slip  boundaries  that  suppressed  the 
formation  of  a  boundary  layer.  The  first  frame  at  0  /zs  shows  the  initial  conditions,  a  planar  incident  shock 
I  moving  from  the  left  to  the  right,  and  the  flame  F.  The  sequence  of  frames,  with  time  increasing  from  top 
to  bottom,  shows  the  incident  shock  I  moving  through  and  distorting  the  flame,  reflecting  from  the  back 
wall,  and  then  the  reflected  shock  R2  passing  through  and  further  distorting  the  flame.  Hot  spots  form  in 
unreacted  material,  and  eventually  two  of  the  hot  spots  undergo  transitions  to  detonations  at  593  /is  and 
608  /zs.  The  evolving  structure  of  this  flow  has  many  similarities  to  the  acetylene-air  simulations  shown  in 
Figure  2.1. 

The  sequence  on  the  right  shows  the  development  of  the  flow  for  the  same  mixture  and  value  of  Ms ,  but 
now  for  the  case  with  a  no-slip  wall  that  allows  a  boundary  layer  to  form.  The  boundary  layer  behind  the 
incident  shock  is  about  0.1  cm  thick  and  is  not  visible  in  these  frames.  From  the  beginning  of  the  simulation 
to  about  225  /zs,  the  flow  development  is  practically  the  same  as  it  was  for  the  case  with  free-slip  walls. 
Boundary-layer  effects  do  become  apparent  at  171  /zs  when  R1  begins  to  form  a  weak  bifurcated  structure 
Bl,  but  Bl  does  not  affect  the  flame  until  a  later  time.  Another  bifurcated  structure  B2  (not  shown)  begins  to 
grow  immediately  after  the  incident  shock  reflects  from  the  back  wall,  but  it  is  quickly  destroyed  by  a  strong 
oblique  shock  coming  out  of  the  flame.  Multiple  reflections  of  oblique  shocks  and  rarefaction  waves  inside 
the  thin  layer  of  unreacted  material  between  the  flame  surface  and  the  bottom  wall  prevent  the  formation 
of  new  bifurcated  structures  until  about  340  /zs. 

The  flame  remains  essentially  unaffected  by  boundary-layer  phenomena  until  363  /zs.  By  that  time, 
R2  has  passed  through  the  flame,  and  a  new  strong  bifurcated  structure  B3  begins  to  develop.  When  a 
bifurcated  foot  forms  and  grows,  a  recirculation  region  forms  and  grows  behind  it.  Near  the  bottom  wall,  the 
flow  in  the  recirculation  region  is  moving  towards  the  shock.  The  flame,  which  almost  touches  the  bottom 
wall  and  moves  with  the  flow  behind  the  shock,  penetrates  the  recirculation  region.  The  flame  approaches 
the  shock  first  as  a  thin  jet  along  the  wall,  and  then  quickly  spreads  inside  the  recirculation  zone.  The  energy 
released  by  the  flame  in  the  recirculation  region  accelerates  the  growth  of  the  bifurcated  foot  until  the  top 
of  the  foot  reaches  the  symmetry  plane  (top  boundary  of  the  computational  domain)  and  a  forms  a  Mach 
reflection  at  the  top  boundary  by  424  /zs. 

The  surface  area  of  the  flame  interacting  with  the  bifurcated  shock  increases,  thus  increasing  the  energy- 
release  rate  in  the  system.  The  energy  release  gradually  increases  the  strength  of  the  bifurcated  shock  and 
the  Mach  stem.  Eventually,  the  temperature  behind  the  Mach  stem  becomes  high  enough  to  ignite  the 
mixture.  Ignition  produces  several  hot  spots,  spontaneous  waves,  and  residual  flames  that  grow  quickly, 
interact  with  vortices,  and  eventually  trigger  a  detonation  Dl. 
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Figure  2.3.  Caption  on  next  page 
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Figure  2.3.  Three-dimensional  simulations  [46]  of  the  same  problem  shown  in  2D  in  Figure  2.1.  The  figure 
shows  shock-flame  interactions  and  DDT  in  a  stoichiometric  ethylene-air  mixture  at  Po  =  1.44  x  103  Pa 
calculated  for  the  initial  and  boundary  conditions  schematically  shown  on  the  top.  Gray  box  on  the  top 
corresponds  to  the  computational  domain.  No-slip  walls  are  y  =  0  and  z  =  0,  opposite  walls  are  symmetry 
planes.  Frames  show  a  time  sequence  of  computed  3D  surfaces  of  flame  (blue  side  faces  unburned  mixture, 
purple  side  faces  burned  gas)  and  shock  (gray  semi-transparent)  in  a  10.6  cm  long  part  of  the  computational 
domain  adjacent  to  the  end  wall  ( x  =  30.4  cm)  for  217-386  //s,  and  shifted  1.5-3  cm  from  the  end  wall  for 
other  times.  Time  in  microseconds  is  given  in  each  frame.  The  letters  indicate  the  incident  shock  (I),  the 
flame  surface  (F),  the  reflected  shock  (R),  the  bifurcated  shocks  (Bl,  B2),  the  recirculation  areas  (Cl, 
C2),  the  shock- flame  attachment  line  (A),  the  Mach  stems  (Ml,  M2),  and  the  detonations  (Dl,  D2). 

The  incident,  planar  shock  moves  along  the  x— axis  and  interacts  with  the  flame.  As  in  the  two- 
dimensional  calculations  described  above  (Figure  2.2),  the  interaction  is  dominated  by  the  Richtmyer- 
Meshkov  instability,  triggered  by  the  flame  curvature.  The  instability  generates  a  large  three-dimensional 
vortex  that  distorts  the  flame  and  creates  a  funnel  of  unreacted  material  in  central  parts  of  the  shock  tube 
near  symmetry  planes.  On  smaller  scales,  the  RM  instability  creates  small  ripples  on  the  flame  surface.  The 
interaction  also  distorts  the  initially  planar  shock  that  propagates  through  burned  and  unburned  materials 
at  different  speeds.  Oblique  parts  of  the  distorted  shock  reflect  from  side  walls  and  interact  again  with  the 
flame,  creating  more  distortions. 

As  the  flame  continues  to  move  with  the  flow,  the  incident  shock  reaches  the  end  wall  at  255  fis  and 
reflects.  The  stronger  reflected  shock  R  moves  back  through  the  already  perturbed  flame  and  creates  more 
distortions  at  different  scales.  The  flame  quickly  develops  a  very  convoluted  surface  and  becomes  turbulent. 
The  energy  release  accelerates  due  to  the  increased  surface  area  of  the  flame  and  the  increased  temperature 
and  density  of  unburned  material  behind  the  reflected  shock. 

When  the  reflected  shock  crosses  the  flame  and  enters  the  unburned  material  upstream,  it  begins  to 
interact  with  boundary  layers  that  developed  near  the  no-slip  walls.  These  interactions  result  in  shock 
bifurcations  Bl  and  B2  that  appear  as  oblique  shocks  propagating  along  the  walls,  y  =  0  and  2  =  0.  Shock 
bifurcations  also  create  strong  recirculation  regions  Cl  and  C2  near  the  walls.  The  recirculation  flow  entrains 
the  flame  that  spreads  through  the  recirculation  region,  approaches  the  bifurcated  shock,  and  then  attaches 
to  it  near  the  wall.  The  attachments  first  occurs  near  symmetry  planes,  where  the  initial  flame  surface  was 
closest  to  no-slip  walls.  Then  the  attachment  line  approaches  the  corner  y  =  z  =  0  between  these  walls. 
The  attachment  line  is  rather  irregular  because  the  flame  and  the  boundary  layer  behind  the  bifurcated 
shock  are  turbulent.  Above  the  attachment  line,  the  flame  surface  F  coincides  with  a  slip  surface  that  forms 
behind  the  bifurcated  shock.  The  Kelvin-Helmholtz  instability  that  develops  along  this  surface  generates 
more  turbulence. 

The  reactive  bifurcated  shocks  propagate  with  a  speed  of  about  a  half  Dqj  and  interact  with  each  other 
along  the  line  that  begins  in  the  corner  y  =  2  =  0.  The  interaction  between  oblique  shocks  Bl  and  B2 
produces  an  oblique  Mach  stem  Ml  that  generates  additional  vorticity  and  locally  increases  the  temperature 
of  shock-compressed  material.  The  oblique  shocks  are  relatively  weak,  and  therefore  the  oblique  Mach  stem 
is  not  strong  enough  to  ignite  the  material. 

DDT  producing  detonation  D2  occurs  when  the  bifurcated  shock  Bl  attached  to  the  wall  2  =  0  forms 
another  Mach  stem  M2  near  the  symmetry  plane  2  =  1.9,  and  this  central  Mach  stem  interacts  with  the 
oblique  Mach  stem  Ml  creating  a  hot  spot.  Another  detonation  Dl  appeared  earlier  in  the  funnel,  near  the 
end  wall,  but  it  did  not  affect  the  formation  of  D2  behind  two  interacting  Mach  stems. 
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Figure  2.4.  Accelerating  flame  (left  column),  DDT,  and  quasi-detonation  (two  right  columns)  in  a  half¬ 
channel  with  obstacles  computed  for  the  conditions  shown  by  the  schematic  on  top  [7].  d/2  =  2  cm,  the 
total  channel  length  L  =  64  cm,  and  dxmin  =  1/512  cm.  Times  in  milliseconds  are  shown  in  frame  corners. 
The  channel  is  filled  with  a  stoichiometric  hydrogen-air  mixture  at  Po  =  1  atm. 

A  laminar  flame,  ignited  at  the  left  top  corner  of  the  computational  domain,  propagates  with  the  velocity 
close  to  Si  =  298  cm/s  relative  to  unburned  material.  Hot  reaction  products  expand  and  push  unreacted 
material  towards  the  open  end  of  the  channel.  The  flame  front  propagates  with  the  moving  flow  and  quickly 
becomes  very  convoluted  as  the  flow  interacts  with  obstacles.  The  increasing  flame  surface  area  results  in 
faster  energy  release,  thus  accelerating  the  flow  and  increasing  the  flame  speed  in  the  laboratory  frame  of 
reference. 

Caption  continued  on  next  page. 
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Figure  2.4.  Caption  continued 

As  the  flame  passes  obstacles,  it  wrinkles  due  to  the  Rayleigh-Taylor  (RT)  instability  caused  by  the  flow 
acceleration.  The  unreacted  flow  ahead  of  the  flame  becomes  sonic  by  1.4  ms,  just  past  O5.  Noticeable  shocks 
begin  to  form  ahead  of  the  flame  past  O7  at  1.85  ms.  They  reflect  from  obstacles  and  side  walls,  and  interact 
with  the  flame  triggering  Richtmyer-Meshkov  (RM)  instabilities.  Kelvin-Helmholtz  (KH)  instabilities  develop 
at  the  flame  surface  when  a  jet  of  hot  burned  material  passes  through  a  narrow  part  of  the  channel  and  a 
shear  layer  forms  downstream  of  the  obstacle.  RT,  RM,  and  KH  instabilities,  and  flame-vortex  interactions 
in  obstacle  wakes  are  the  mechanisms  responsible  for  increases  in  flame  surface  area,  energy-release  rate, 
and,  eventually,  shock  strength.  The  elevated  temperature  behind  shocks  also  contributes  to  the  increased 
energy-release  rate  because  Si  increases  and  shocks  passing  through  the  reaction  zone  release  additional 
energy. 

The  average  flame  velocity  gradually  increases  and  reaches  800  m/s  by  2.f  ms.  This  velocity  equals  about 
0.8  of  sound  speed  in  the  burned  material,  or  0.4Dcj,  and  is  typical  of  what  has  been  called  the  choking 
regime  of  flame  propagation  observed  in  experiments  with  obstructed  channels.  As  the  shock  and  the  flame 
accelerate,  the  leading  edge  of  the  flame  remains  about  1  cm  behind  the  leading  shock,  which  diffracts  at 
every  obstacle  and  reflects  from  the  bottom  wall  after  each  diffraction.  The  reflection  type  changes  from 
regular  to  strong  as  the  reflection  point  approaches  the  next  obstacle.  The  resulting  Mach  stem  becomes 
stronger  after  each  diffraction,  and  the  temperature  of  the  hot  region  that  forms  when  the  Mach  stem  collides 
with  an  obstacle  increases.  At  2T00  ms,  the  reflection  of  the  Mach  stem  from  O12  creates  a  region  with 
temperatures  above  830  K.  Two  hot  spots  in  this  region  ignite  producing  two  small  flame  kernels.  Then  a 
detonation  appears  near  one  of  the  kernels  and  propagates  through  the  unreacted  material. 

The  newly  formed  detonation  propagates  through  the  gap  between  the  flame  and  the  obstacle  into 
the  shock-compressed  material  ahead  of  the  flame.  As  the  detonation  passes  around  the  obstacle,  the 
lower  part  of  the  front  decouples  into  a  separated  shock  and  a  flame.  The  upper  part  of  the  front  remains 
essentially  undisturbed  and  develops  detonation  cells  before  it  collides  with  the  upper  boundary.  The  collision 
creates  a  strong  reflected  shock  that  triggers  a  detonation  in  both  the  shock-compressed  layer  between  the 
leading  shock  and  the  decoupled  flame,  and  the  uncompressed  material.  The  strong  detonation  wave  in  the 
uncompressed  material  quickly  develops  detonation  cells,  collides  with  O13  at  2.125  ms,  and  diffracts.  As  the 
diffraction  weakens  the  detonation  wave,  detonation  cells  grow  and  form  an  irregular  two-level  structure.  The 
diffraction  on  O14  completely  decouples  the  shock  and  flame  by  2.164  ms,  and  effectively  kills  the  detonation. 
A  new  detonation  is  ignited  in  the  shock-compressed  material  by  the  collision  of  the  Mach  stem  with  O1.5  at 
2.179  ms,  but  this  detonation  is  unable  to  propagate  through  the  very  narrow  gap  between  the  obstacle  and 
the  flame.  The  leading  shock  and  the  flame  remain  decoupled  until  the  Mach  stem  hits  Oi6  and  triggers  a 
new  detonation  at  2.217  ms  that  spreads  past  the  obstacle.  This  quasi-detonation  regime  that  involves  the 
detonation  diffraction,  failure,  and  reignition  continues  until  the  end  of  the  channel. 
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Figure  2.5.  Flame  and  shock  configurations  in  hydrogen-air  mixture  just  before  the  detonation  initiation 
or  at  the  end  of  the  channel.  A  Mach  stem  forms  near  the  bottom  wall  for  larger  channels  ( S  =  4,  6,  8 
cm),  but  not  for  smaller  channels  ( S  =  1,  1.5,  2,  3  cm).  Half-channel  height  d/2=2  cm.  S  is  defined  by  the 
schematic  in  Figure  2.4. 
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(a)  Initial  conditions 


(b)  (4000)  0.21  ms 


(c)  (30,000)  1 .59  ms 


Figure  2.6  Premixed  hydrogen-air  flame  propagating  through  2D  array  of  obstacles  after  a  weak  spark 
ignition  in  the  left  bottom  corner.  The  size  of  each  obstacle  is  1  x  1  cm2,  the  distance  between  neighboring 
rows  of  obstacles  is  1  cm,  and  the  computational  cell  size  1/128  cm.  The  timestep  number  and  time  from 
ignition  are  shown  above  each  frame.  Colors  show  temperature.  Detonations  appear  in  frame  (h)  and  spread 
in  all  directions.  The  dark  (low-temperature)  band  in  frame  (i)  separates  materials  burned  in  deflagration 
and  detonation  regimes.  An  enlargement  of  a  portion  of  frame  (g)  was  shown  earlier  in  Figure  1.2. 
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3.  Flame  Acceleration,  DDT,  and  Detonations  in  Methane- Air  Mixtures 

The  purpose  of  the  theoretical  and  numerical  effort  is  to  develop  a  tool  that  will  help  us  to  address  the 
pivotal  question: 

Given  a  flammable  mixture  of  methane  and  air  of  sufficient  volume  and  lateral  extent,  can  the  mixtures 
develop  into  a  detonation  from  a  weak  spark  ignition  source? 

To  address  this  question,  we  defined  four  tasks  that,  in  conjunction  with  well-defined  experiments,  will  allow 
us  to  answer  the  question  and  give  us  the  basis  for  understanding  the  results  as  well  as  extrapolating  them 
to  new  situations.  In  summary,  these  tasks  are: 

1.  Develop  and  test  a  model  for  stoichiometric  methane- air  combustion.  Perform  a  series  of  simulations 
of  flame  propagation,  acceleration,  and  possible  DDT  in  channels  with  obstacles. 

2.  Develop  and  test  similar  models  for  fuel-rich  and  fuel-lean  combustion. 

3.  Examine  the  effects  of  spatially  varying  fuel  concentrations  in  a  simulation. 

4.  Document  the  effects  of  stochasticity,  i.e. ,  of  large-scale  uncertainty  in  system  behavior  caused  by 
minute  perturbations  in  initial  conditions  or  numerical  procedures,  on  flame  acceleration  and  DDT. 

The  progress  in  each  of  these  areas  is  described  below.  In  this  exposition,  we  have  moved  the  discussion  of 
Task  4  to  just  after  Task  1  because  the  tests  were  done  on  the  stoichiometric  mixture. 

The  final  model  developed  must  have  the  ability  to  compute  the  different  stages  of  the  evolution  of  a 
chemically  reactive  flow:  ignition  by  a  small  spark,  rapid  flame  acceleration,  development  of  shocks,  shock- 
flame  interactions,  and  detonation  initiation  in  complex  geometries  for  natural  gas  mixtures  with  spatially 
and  temporally  varying  stoichiometries.  The  solutions  must  come  from  the  set  of  coupled,  reactive  Navier- 
Stokes  equations  (Eqs.  2.1  -2.5).  What  we  now  describe  are  the  first  steps  toward  achieving  this  and  a  reliable 
predictive  capability  for  DDT. 

3.1.  DDT  in  Stoichiometric  Methane- Air  Mixtures  (Task  1) 

The  items  in  Task  1  are  the  first,  necessary  steps  in  developing  a  multidimensional  numerical  model.  The 
material  presented  here  is  described  in  more  detail  in  references  [9-13]. 

3.1.1.  Development  of  the  Chemical  Model:  Parameter  Calibration 

Numerical  models  that  can  describe  the  behavior  of  shocks  and  detonations  vary  widely  in  their  complex¬ 
ity,  but  for  many  practical  situations,  an  extensive  description  of  the  details  of  the  chemical  pathways  is 
expensive  and  unnecessary.  Instead,  it  is  more  important  to  have  an  accurate  model  of  the  fluid  dynamics 
coupled  to  a  model  for  the  chemical-energy  release  that  puts  the  released  energy  in  the  “right”  place  in 
the  flow  at  the  “right”  time.  For  detonations,  direct  shock  ignition  and  subsequent  propagation,  as  well  as 
other  more  complex  transient  detonation  structures,  can  be  “accurately  enough”  computed  if  an  acceptable 
representation  of  the  chemical  induction  time  (i.e.,  autoignition  time)  is  known  as  a  function  of  the  state 
variables,  temperature,  and  pressure  [50].  Flame  propagation  and  structure  can  be  “accurately  enough” 
computed  with  a  model  of  energy  release  and  diffusion  coefficients  that  create  a  reasonable  preheat  zone 

[51]- 
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The  issue  is  always:  “What  is  enough?”  How  accurate  or  detailed  does  the  model  have  to  be?  This 
question  has  led  to  a  series  of  simplified,  single-step  models  for  either  flames  or  detonations  [51  -54]  separately. 
What  we  need  for  prediction  of  DDT,  however,  is  more  demanding.  We  need  a  single  model  that  will  work 
for  both  combustion  regimes.  In  previous  work,  we  described  similar  single-step  models  for  low-pressure 
acetylene,  low-pressure  ethylene,  and  atmospheric-pressure  hydrogen-air  mixtures  [5-7,  41,  43,  44,  55].  Here 
we  address  the  problem  of  developing  a  minimal  model  that  captures  the  essential  features  of  a  flame,  a 
detonation,  and  the  deflagration-to-detonation  transition  (DDT)  in  nrethane-air  mixtures. 

Consider  the  physical  phenomena  we  need  to  model.  Flame  acceleration  in  the  laboratory  frame  occurs 
primarily  due  to  advection  by  the  induced  gas  flow,  which  can  be  orders  of  magnitude  larger  than  the  laminar 
and  turbulent  flame  speeds  [6,  7].  The  gas  flow  is  driven  by  thermal  expansion  of  the  combustion  products 
and  increases  with  the  amount  of  heat  released  in  a  unit  of  time  by  the  flame  front,  which  is  a  function  of  the 
surface  area  of  the  flame.  The  flame  surface  area  increases  due  to  stretching  by  the  flow  and  wrinkling  caused 
by  turbulent  motions  and  fluid-dynamic  instabilities.  When  the  flow  speed  approaches  the  speed  of  sound, 
shocks  form,  and  shock-flame  interactions  become  an  important  mechanism  for  the  flame  wrinkling  and 
turbulence  generation.  Hot  spots,  which  can  evolve  to  generate  flames,  shocks,  and  detonations,  may  arise 
in  unreacted  material  from  many  types  of  interactions  in  a  turbulent,  shock-laden  reactive  flow.  In  these  hot 
regions,  a  spatial  temperature  gradient  exists,  inside  of  which  the  temperature  can  be  high  enough  to  ignite 
the  reactants.  Ignition  occurs  consecutively  in  multiple  layers  of  material  heated  to  different  temperatures, 
thus  forming  a  reaction  wave  [40] .  This  wave  can  generate  a  strong  shock  and  eventually  a  detonation.  The 
survival  of  the  newly  formed  detonation  wave  then  depends  on  local  thermal  and  chemical  conditions  and 
geometrical  constraints  [5]. 

Based  on  this  progression  of  physical  processes,  the  properties  of  the  fuel-air  mixture  that  are  important 
during  the  flame-acceleration  stage  are  the  laminar  flame  speed  Si,  the  adiabatic  flame  temperature  T/,,  the 
viscosity  v,  the  laminar  flame  thickness  Xi,  and  the  speed  of  sound  c.  In  the  work  that  follows,  we  assume 
that  the  gas  mixtures  behave  as  ideal  gases,  which  means  pressure  P  and  internal  energy  e  are  related  to 
the  density  and  temperature  by  the  state  equations 


p=p^L 

M 


3.1 


P 

e  (t^T)  • 


3.2 


Here  the  internal  energy  is  e  =  E  —  pv2/ 2,  where  E  is  the  total  energy  and  v  is  the  velocity.  The  sound 
speed  is  given  by  c  =  RT,  and  the  specific  heat  capacity  at  constant  pressure  is  related  to  the  gas 
constant  by  cp  =  7??/(7  —  1).  The  temperature  rise  in  an  adiabatic  system  due  to  chemical-energy  release 
is  Tb  —  Ta  =  q/cp  =  q( 7  —  l)/7 R.  Thus,  the  adiabatic  flame  temperature  depends  on  the  initial  system 
temperature  T0,  the  heat  release  q,  and  the  ratio  of  specific  heats  7.  The  laminar  flame  speed  also  depends 
on  q  and  7,  but  it  is  also  strongly  influenced  by  the  rate  of  chemical  energy  release,  which  is  determined  by 
the  choice  of  reaction  model.  Here  we  assume  that  fuel,  F,  is  converted  to  product,  P,  via  a  single  reaction, 
F  — >  P,  that  has  an  Arrhenius  dependence  on  temperature,  characterized  by  a  pre-exponential  factor  A  and 
an  activation  energy  Ea.  The  thermal  diffusivity  of  the  reactants,  k  =  K/pcp,  determines  the  length  of  the 
preheat  zone  of  the  laminar  flame,  and,  thus,  affects  both  Si  and  Xi. 
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The  late  stages  of  DDT  are  controlled  by  induction  delays  behind  strong  shocks  that  correlate  with 
some  properties  of  detonation  waves.  Detonations  are  inherently  unstable,  leading  to  the  formation  of  highly 
localized  regions  of  high  pressure  and  vorticity.  Three  shocks  intersect  to  form  the  backbone  of  these  regions, 
which  are  commonly  referred  to  as  triple  points.  Along  the  leading  edge  of  the  reaction  front,  one  shock 
is  a  detonation  wave  and  the  second  is  an  inert  shock  that  has  decoupled  from  the  reaction  front.  A  third 
transverse  shock  is  inclined  to  the  two  leading  shocks  behind  the  leading  edge.  That  shock  can  either  be 
inert  or  reactive  (a  detonation).  Several  examples  of  triple  points  are  shown  below  in  Figure  3.5.  The  triple 
points  move  along  the  detonation  front  in  the  direction  transverse  to  the  direction  of  detonation  propagation. 
The  number  and  trajectory  of  triple  points  are  characteristic  of  the  properties  of  the  gas  mixture.  When 
multiple  triple  points  form  along  a  detonation  front,  they  leave  behind  a  cellular  pattern  that  has  been 
observed  in  experiments.  Since  induction  delays  correlate  with  the  size  of  these  cellular  structures,  they  can 
be  used  to  calibrate  a  model  that  should  approximately  describe  key  phenomena  responsible  for  the  final 
stages  of  DDT.  The  detonation  properties  used  for  the  calibration  are  the  theoretical  Chapman-Jouget  (CJ) 
detonation  velocity  Dqj  and  the  detonation  cell  size  A.  The  velocity  Dqj  depends  only  on  the  total  amount 
of  heat  release  and,  hence,  on  q  and  7.  If  the  sound  speed  is  correct,  Dqj  gives  an  indication  of  the  correct 
value  of  q ,  as  X),  does  for  a  deflagration.  Ignition  delay  times  also  depend  on  the  rate  of  chemical  energy 
release.  These  delay  times  correlate  well  with  the  reaction  zone  structure,  particularly  with  the  characteristic 
size  of  the  detonation  cells,  A,  which  is  controlled  by  q,  7,  A  and  Ea. 

The  form  of  the  model  we  are  seeking  for  the  reaction  rate  Cl  is 

Cl  =  ApYexp{—Ea/RT),  (3.3) 

where  Y  is  the  fuel  mass  fraction.  The  transport  coefficients,  viscosity  v,  mass  diffusivity  D ,  and  thermal 
diffusivity  k  vary  with  temperature  according  to 

rj-in  rj-ifl  rpn 

V  =  v0 — ,  D  =  Do  — ,  K  =  Ko  —  •  (3.4) 

p  p  p 

The  parameters  17,,  D0,  and  n0  are  assumed  to  be  constant,  and  n  is  a  constant  chosen  to  be  0.7.  In  this 
model,  the  ratio  of  specific  heats  7  does  not  vary  with  temperature.  Target  output  for  the  quantities  Si,  Tb , 
Xi,  Dcj ,  Xd,  and  A  for  a  stoichiometric  (9.5%)  nrethane-air  mixture  are  given  in  Table  1. 


Table  1.  Target  Values  for  Output  of  Methane- Air  Models 


Quantity 

Lean 

(7%,  0  =  0.72) 

Stoichiometric 

(9.5%,  0=1) 

Rich 

(13%,  0  =  1.43) 

Si 

15-24  cm/s 

34-45  cm/s 

15-25  cm/s 

Tb 

1850-1940  K 

2200-2330  K 

1960-2000  K 

Xl 

0.17-0.2  cm 

0.08-0.12  cm 

0.18-0.25  cm 

Dcj 

1680  m/s 

1815  m/s 

1800  m/s 

Xd 

0.8-1. 6  cm 

0.13-0.62  cm 

0.69-1.1  cm 

A 

«  80  cm 

13-31  cm 

ss  45  cm 
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The  model  calibration  is  started  by  performing  a  series  of  one-dimensional  calculations  of  flame  and 
detonation  structures  for  a  range  of  input  parameters.  Then  we  choose  a  set  of  parameters  that  most  closely 
reproduces  both  the  laminar  flame  properties  and  the  detonation  properties.  After  this,  two-dimensional 
simulations  are  used  to  compute  detonation  cell  sizes  using  the  model  parameters  arrived  upon  in  the  previous 
step.  Finally,  full  simulations  using  the  model  parameters  are  used  to  compute  flame  acceleration  and  DDT, 
and  these  results  are  compared  to  experiments. 

First,  the  values  of  t'o  and  a'o  were  calculated  directly  from  the  known  viscosity  and  thermal  conductivity 
of  air  at  298  K  and  1  atm.  For  simplicity,  we  assumed  that  the  Lewis  number  (Le  =  Kq/Dq)  of  the  mixture 
is  equal  to  unity,  i.e. ,  that  mass  and  heat  diffuse  at  the  same  rate,  so  that  Dq  =  no-  This  approximation  is 
appropriate  for  mixtures  of  methane  and  air  for  which  k0  is  only  slightly  larger  than  Da.  Assuming  Le  =  1 
greatly  simplifies  the  one-dimensional  analyses  used  to  perform  the  model  calibration.  We  then  computed 
properties  of  a  steady,  one-dimesional  laminar  flame  by  solving  an  ordinary  differential  equation  describing 
thermal  conduction  and  energy  release  inside  a  steady-state  reaction  wave, 

f  =  ,(t/Cpf-»o)  (3.5) 


Ft  =  K —  (3.6) 

ax 

where  U  =  Sip^/p  is  the  flow  velocity  in  the  frame  moving  with  the  flame,  Ft  is  the  thermal  flux,  and  Q  is 
the  reaction  rate  given  by  Eq.  (3.3).  According  to  Eqn.  (3.4),  the  thermal  conductivity  K  is  a  function  of 
temperature,  K  =  npcp  =  kqT0'7 pcp.  Details  of  the  solution  method  are  given  in  [9,  56].  We  search  for  a  set 
of  parameters,  A,  q,  Ea,  and  7,  for  which  the  computed  value  of  Si  matches  experimental  data  [57,  58],  and 
that  of  Tj,  matches  values  calculated  using  a  complex  reaction  mechanism  for  methane-air  combustion  [33]. 
The  laminar  flame  thickness,  xi,  is  based  on  the  temperature  gradient,  i.e.,  Xi  =  (TJ,  —  To)/ max \dT/dx\. 

We  use  a  Zeldovich-von  Neumann-Doering  (ZND)  model  to  compute  the  half-reaction  thickness  of  one¬ 
dimensional  detonations  using  the  same  reaction  model  described  above  (Eq.  3.3)  for  a  set  of  parameters  A, 
Ea,  7,  and  q.  The  reaction  zone  of  a  one-dimensional  detonation  is  described  by 


dp  _  q  ft  p  (7  -  1) 
dt  (U*2  —  c2) 

(3.7) 

de  P  dp  ■ 

dt  p2  dt  ^ 

(3.8) 

dx* 

-  =  U 

dt 

(3.9) 

where  e  is  the  internal  energy,  t  is  time,  x*  is  the  distance  from  the  shock,  U* 

=  DcjPo/p  is  the  flow  velocity 

in  the  shock  frame,  c  =  y/'yP/p  is  the  sound  speed,  and 


F>cj 


Co 


g  poh2  - 1) 
p0  27 


q  p0(72  -  l)\ 

Po  27  J 


(3.10) 


is  the  Chapman-Jouget  detontation  velocity.  The  solution  procedure  for  this  set  of  equations  is  described  in 
detail  in  [9  56]. 
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Solutions  of  Eqns.  (5)-(8)  give  the  reaction-zone  profile,  from  which  we  find  the  half-reaction  thickness 
of  the  one-dimensional  detonation  wave,  Xd-  This  quantity  correlates  with  the  detonation  cell  size  A  [59], 
which  is  often  measured  in  detonation  experiments.  Even  though  detonation  cells  do  not  develop  until  after 
DDT  occurs,  A  is  generally  used  instead  of  Xd  in  empirical  correlations  related  to  the  detonation  initiation 
and  DDT  [60-71]  For  example,  experimental  evidence  suggests  that  the  ratio  of  the  system  size  to  A  must 
be  greater  than  one  for  a  sustained  detonation  to  occur  [60].  We  use  measured  values  of  A  [72]  (and, 
equivalently,  Xd)  and  values  of  Dqj  calculated  using  a  thermodynamic  equilibrium  code  (e.g.,  [73])  to  find 
a  set  of  consistent  model  parameters. 

The  choice  of  parameters  7  and  q  determines  Tb  for  laminar  flames  and  Dqj  for  one-dimensional 
detonations.  Curves  representing  values  of  q  that  give  Tb  =  2210  K  and  Dqj  =  1820  m/s  as  a  function 
of  7  are  shown  in  Figure  3.1a.  We  choose  the  values  of  q  and  7  at  the  intersection  of  these  two  curves  to 
use  in  our  model.  The  remaining  parameters,  Ea  and  A,  determine  Si  and  Xd  for  fixed  ko,  7  and  q.  The 
determination  of  A  for  such  a  high-activation-energy  mixture  is  not  precise  because  the  detonation  cells  are 
highly  irregular.  Hence,  ranges  for  A  and  Xd  are  given  in  Table  1.  A  representative  value  within  the  range 
shown  in  the  table,  Xd  =  0.229  cm,  was  used  as  the  target  value  for  the  calibration  process.  Figure  3.1b 
shows  the  values  of  A  for  which  Si  =  38.02  cm/s  and  Xd  =  0.229  cm  as  a  function  of  Ea.  Again,  we  choose 
A  and  Ea  at  the  point  of  intersection  of  these  curves.  Composite  model  values  with  their  output  quantities 
are  shown  in  Table  2.  We  note,  however,  that  a  range  of  Ea  and  A  exists  that  would  give  satisfactory  values 
of  A. 


Table  2.  Material,  chemical,  and  reaction-wave  parameters  for  stoichiometric  (9.5%)  methane- 
air 


Quantity 

Value 

Definition 

Input 

Po 

1  atm 

Initial  pressure 

T0 

298  IC 

Initial  temperature 

Po 

1.1042  x  10"3  g/cm3 

Initial  density 

7 

1.197 

Adiabatic  index 

Ms 

27 

Molecular  weight 

A 

1.64  x  1013  cm3/g/s 

Pre-exponential  factor 

Ea 

67.55  R  To 

Activation  energy 

q 

39.0  R  To  /M 

Chemical  energy  release 

vo 

3.6  x  10-6  (g/s-cm-K0-7) 

Transport  constants 

0 

II 

to 

0 

6.25  x  10"6  (g/s-cm-K0-7) 

Transport  constants 

Output 

Si 

38.02  cm/s 

Laminar  flame  speed 

Tb 

2210  K 

Adiabatic  flame  temperature 

Xl 

0.0439  cm 

Laminar  flame  thickness 

Dcj 

1820  m/s 

Chapman-Jouguet  detonation  velocity 

Xd 

0.229  cm 

ID  half-reaction  detonation  thickness 
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3.1.2.  Description  of  Numerical  Method  and  Validation  of  Basic  Model 

In  order  to  properly  compute  DDT,  a  numerical  model  must  be  able  to  compute  a  wide  variety  of 
combustion  waves,  ranging  from  low-speed  laminar  flames  to  moderate-speed  turbulent  flames  to  high-speed 
detonations.  In  this  section,  we  explore  how  the  model  developed  in  the  previous  section  can  be  used  to 
compute  the  limiting  cases  of  this  range  of  combustion  waves:  a  laminar  flame  and  a  detonation.  Here  we 
present  computations  of  two-dimensional  unsteady  laminar  flames  and  detonations  and  compare  them  to  the 
idealized,  steady-state  ID  flame  and  detonation  structures  computed  in  Section  3.1.1. 

The  code  used  to  perform  these  validation  tests  and  all  subsequent  simulations  described  in  Sections 
3  and  4  solves  the  RNSE  using  an  explicit,  second-order,  Godunov-type  numerical  scheme  incorporating  a 
Riemann  solver.  The  domain  is  discretized  using  a  structured  adaptive  mesh  based  on  the  fully  threaded 
tree  data  structure  [74].  The  mesh  refinement  is  dynamically  controlled  by  gradients  of  density,  temperature, 
and  composition. 

3.1.2a.  Laminar  Flame  Structure 

Equations  2.1  -2.5  in  conjunction  with  the  calibrated  reaction  model  (Eqn.  3.3  and  Table  2)  were  solved 
using  the  standard  explicit,  second-order,  numerical  scheme  described  above.  We  consider  a  thin  (in  the 
y-direction),  two-dimensional  channel  that  is  initially  filled  with  premixed  methane  and  air.  An  initial 
approximation  for  a  planar  laminar  flame  is  placed  across  the  domain  so  that  the  reaction  wave  propagates 
into  the  fresh  reactants  in  the  positive  x-direction.  After  a  short  adjustment  period,  the  correct  laminar 
flame  profile  forms,  and  the  flame  propagates  at  a  constant  velocity,  Si.  Here  we  compare  the  computed  flame 
structure  and  propagation  speed  to  those  obtained  with  the  ID  steady-state  approximation  described  earlier. 
We  consider  two  different  grid  resolutions  at  the  flame  front,  dxmin  =  0.018125  cm  («  3  computational  cells 
per  laminar  flame  thickness)  and  dxmin  =  0.00113  cm  («  39  cells/x;).  Away  from  shocks  and  flame  fronts, 
the  maximum  computational  cell  size  was  taken  to  be  dxmax  =  0.29  cm. 

Figure  3.2  shows  the  ID  flame  structure  calculated  for  these  two  different  grid  resolutions  along  with  the 
ID  steady-state  solution  (lines).  The  high-resolution  calculation  reproduces  the  theoretical  flame  structure 
and  gives  nearly  the  same  flame  speed,  Si  =  38.02  cm.  There  are  some  differences  in  the  reaction-rate  profiles 
calculated  at  lower  resolution  ( dxmin  =  0.018125  cm).  The  computed  laminar  flame  speed  is  approximately 
12%  smaller  than  the  theoretical  laminar  flame  speed.  More  resolution  tests  are  discussed  below  in  Task 
4.  Here  we  only  note  that  when  we  compute  DDT  in  larger  channels,  dxmin  is  limited  by  the  available 
computational  resources,  and  in  most  cases,  the  simulations  must  be  somewhat  under-resolved.  The  coarser 
resolution  simulation  shown  in  Figure  3.2  was  chosen  to  reflect  the  minimum  amount  of  grid  resolution  we 
used  in  the  subsequent  DDT  simulations. 

3.1.2b.  Detonation  Structure 

Multidimensional  instability  of  gas-phase  detonations  results  in  a  complex,  dynamic  structure  at  a 
detonation  front.  This  structure  involves  multiple  transverse  waves  and  triple-shock  configurations  that 
interact  with  each  other,  causing  the  energy  release  to  be  highly  nonuniform  [75].  The  size  and  regularity  of 
a  detonation  cell  characterize  the  kinetics  of  energy  release:  the  cell  size  correlates  with  the  energy-release 
rate  and  the  regularity  is  controlled  by  the  effective  activation  energy  Ea/ RTznd,  where  Tznd  is  the  post¬ 
shock  temperature  in  the  ZND-detonation  structure.  An  example  of  a  very  regular  structure  formed  in  a 
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low-pressure  mixture  highly  diluted  by  an  inert  gas  is  shown  in  [75]. 

Some  reactive  systems,  such  as  metlrane-air,  show  even  more  complex  detonation  cell  structures.  For 
example,  fine  cellular  structures  may  appear  embedded  inside  main  detonation  cells.  In  some  cases,  such  fine 
cellular  structures  may  be  caused  by  two  distinct  stages  of  energy  release  that  occur  at  different  timescales 
[76,  77].  In  many  cases,  however,  multilevel  cellular  structures  can  be  explained  by  secondary  instabilities  of 
overdriven  parts  of  the  detonation  front  that  can  even  appear  in  systems  with  a  one-step  energy  release  and 
high  Ea/ RTznd-  Analysis  of  experimental  data  [78]  and  numerical  simulations  [79-81]  show  that  secondary 
detonation  cells  appear  for  Ea/ RTznd  >  6.5. 

We  have  now  performed  a  series  of  extremely  lrigh-resolution  unsteady  numerical  simulations  of  deto¬ 
nations  propagating  in  a  stoichiometric  methane-air  mixture.  This  is  a  very  high-activation  energy  system, 
Ea/ RTznd  =  17.2,  and,  accordingly,  the  simulations  show  the  formation  of  multiple  levels  of  cellular  struc¬ 
tures  caused  by  the  propagation  of  very  irregular  triple-shock  configurations  [10].  The  simulations  show 
two  main  types  of  triple-shock  configurations:  one  in  which  there  was  a  triple-shock  configuration  with  an 
attached  transverse  detonation,  and  another  triple-shock  configuration  with  an  inert  transverse  wave.  These 
combine  to  produce  three  types  of  collisions,  which  seem  characteristic  of  a  system  of  highly  irregular  struc¬ 
tures  with  very  high-activation  energy.  Here  we  describe  these  simulations  in  more  detail.  Later  in  this 
report  (Section  4)  we  relate  the  results  of  these  simulations  to  recent  experimental  cell  measurements  at 
NIOSH  Lake  Lynn  Laboratory  (LLL). 

The  system  was  modeled  by  solving  the  same  set  of  equations  with  the  same  AMR  approach  to  adaptive 
gridding  as  was  used  to  perform  the  2D  flame  calculations  in  Section  3. 1.2. a  and  all  of  the  DDT  simulations 
described  below.  It  has  been  shown  that  such  models  capture  most  of  the  key  features  responsible  for 
the  formation  of  detonation  cells,  including  the  effects  of  chemical  kinetics  on  cell  regularity  and  formation 
of  secondary  cell  structures  [79-81].  We  consider  a  two-dimensional  rectangular  channel  of  height  32  cm 
and  length  1024  cm.  Using  2D  simulations  allows  us  to  achieve  a  very  high  numerical  resolution  in  the 
current  calculations  that  would  be  prohibitively  expensive  in  3D.  Based  on  previous  experiments  [72]  and 
calculations  [59],  the  channel  height  is  large  enough  to  capture  at  least  one  detonation  cell.  A  significant 
number  of  computational  grid  cells  was  needed  to  resolve  the  instabilities  in  the  detonation  front  properly. 
The  smallest  grid  size  used  in  these  computations  was  1/1024  cm,  which  is  equivalent  to  about  234  cells 
per  half  reaction  thickness.  Details  of  the  ignition  method  is  given  in  [10].  Figure  3.3  shows  the  speed  of 
the  leading  edge  of  the  computed  detonation  front  as  a  function  of  time.  The  average  propagation  speed  of 
1813  m/s  is  consistent  with  the  CJ  detonation  velocity  (1820  m/s),  as  defined  by  thermodynamic  parameters. 

Locally,  the  shape  and  the  velocity  of  the  detonation  front  are  dictated  by  the  dynamics  of  primary 
triple-shock  configurations  that  form  as  a  result  of  the  inherent  instability  of  a  planar  detonation.  At  times, 
the  leading  edge  of  the  front  is  an  overdriven  detonation  with  velocity  greater  than  Dqj,  and  at  other  times 
a  decoupled  shock- flame  complex  with  velocity  less  than  Dqj-  A  sequence  of  density  fields  taken  from  the 
calculations  at  several  instants  over  a  characteristic  detonation-cell  crossing  time  are  shown  in  Figure  3.4. 

The  various  types  of  triple-shock  configurations  appearing  in  the  computation  are  analyzed  in  detail  in 
[10]  and  are  shown  in  Figure  3.5.  Most  of  the  primary  triple-shock  configurations  we  see  in  the  simulations 
involve  transverse  detonations  that  propagate  through  the  large  induction  zone  (labeled  3  in  the  figure) 
behind  weak  parts  of  the  leading  shock  (1).  A  close-up  view  of  the  temperature  field  in  the  vicinity  of  a 
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typical  triple-shock  configuration  is  shown  in  Figure  3.5a.  The  transverse  detonation  wave  (2)  is  nearly 
perpendicular  to  the  leading  shock  (1),  and  burns  most  of  the  unreacted  material  (3)  behind  it.  Even  the 
material  in  the  small  triangular  region  between  the  oblique  shock  (4)  and  the  slip  line  (5)  detonates,  in 
contrast  to  the  case  with  lower  Ea/ RTznd  computed  in  [80],  where  it  did  not  burn  and  left  behind  an 
unreacted  tail. 

Occasionally,  we  saw  a  different  type  of  triple-shock  configuration:  a  transverse  detonation  failed  and 
became  an  inert  transverse  shock  (9).  One  of  these  configurations  is  shown  in  Figure  3.5b.  Here,  the 
material  detonates  only  in  the  overdriven  part  of  the  leading  front  (6).  A  variation  of  this  situation  is  shown 
in  Figure  3.5c  where  only  a  small  part  of  the  transverse  detonation  survives.  A  thin  layer  of  material  is 
burned  in  this  region  in  addition  to  the  burning  in  the  overdriven  part  of  the  leading  front.  These  types  of 
configurations  leave  behind  large  pockets  of  compressed,  unburned  material  (8)  that  may  become  completely 
surrounded  by  burned  gas  when  two  triple  points  collide,  which  can  limit  the  total  amount  of  energy  released 
near  the  reaction  front  and  limit  the  propagation  velocity  of  the  detonation. 

Since  both  types  of  triple-shock  configurations  shown  in  Figure  3.5  exist  in  the  same  system,  three  types 
of  collisions  between  triple-shock  configurations  are  possible.  The  examples  from  the  current  simulation  are 
shown  in  Figure  3.6.  Isolated  unreacted  pockets  form  when  at  least  one  of  the  colliding  transverse  waves  is 
inert.  If  another  transverse  wave  is  a  detonation  (Figure  3.6b),  then  it  continues  to  propagate  into  the  pocket 
after  the  collision  and  burns  the  compressed  material.  If  both  transverse  waves  are  inert  (Figure  3.6a),  then 
the  pocket  drifts  downstream  and  burns  as  a  deflagration.  Collisions  between  two  triple-shock  configurations 
with  reactive  transverse  shocks  (Figure  3.6c)  result  in  significantly  weaker  reflected  triple  shocks  since  there 
is  very  little  unburned  material  behind  the  leading  edges  of  detonation. 

Due  to  the  high  activation  energy  in  the  model  used  in  these  calculations,  both  transverse  detonations 
and  the  overdriven  parts  of  the  leading  front  are  unstable.  They  produce  the  smaller  triple-shock  configura¬ 
tions  that  can  be  seen  in  Figure  3.6.  For  overdriven  parts  of  the  leading  front,  these  secondary  instabilities 
increase  in  size  with  the  distance  from  the  primary  triple  point  (Figure  3.6a).  The  spacing  between  the  sec¬ 
ondary  triple  points  varies  from  fractions  of  a  millimeter  to  several  millimeters.  The  instability  of  transverse 
detonations  occurs  at  even  smaller  scales  of  the  order  of  0.1  mm  or  less  that  are  not  well  resolved  in  the 
simulations. 

Finally,  we  emphasize  that  the  unstable  detonation  front  modeled  here  is  very  irregular.  Even  though 
characteristic  length  scales  associated  with  primary  triple  points  can  be  approximately  estimated,  the  dis¬ 
tinction  between  primary  and  secondary  triple-shock  configurations  often  becomes  blurred.  The  detonation 
instabilities  can  be  observed  for  all  range  of  scales,  from  about  a  half  of  the  channel  width  down  to  the  grid 
resolution.  This  leads  to  the  formation  of  a  very  complicated  detonation  cell  structure.  While  it  is  difficult 
to  make  quantitative  assessments  of  cell  sizes  and  regularity,  the  numerical  results  obtained  using  the  reac¬ 
tion  model  described  in  Section  3.1.1  qualitatively  agree  with  experimental  results.  Direct  comparisons  with 
experiments  performed  at  LLL  will  be  discussed  in  a  later  section. 
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3.1.3.  Channels  with  Obstacles 


We  next  use  Eqns.  (2.1  -2.5)  to  calculate  multidimensional,  turbulent,  accelerating  flames  and  subsequent 
DDT.  Obstructed  channels  promote  faster  flame  acceleration  than  smooth  channels  and,  hence,  provide  a 
useful  geometry  in  which  to  study  DDT.  In  addition,  this  configuration  has  been  used  for  several  experimental 
studies,  and  so  there  is  data  that  can  be  used  for  comparisons.  A  schematic  of  the  channel  used  in  the 
calculations  is  shown  in  Figure  3.7. 

We  consider  three  different  configurations  chosen  to  be  similar  to  experimental  systems  [72,  82]  and 
summarized  in  Table  3.  For  each  test  case,  the  channel  is  uniformly  filled  with  a  stoichiometric  methane-air 
mixture. 

-  Configuration  7.6.  Here  L  =  216.2  cm,  d  =  7.6  cm,  which  models  the  7. 6x7.6  cm  square  channel  used  in 
[82].  In  the  experiments,  both  ends  of  the  channel  were  closed,  and  the  initial  pressure  in  the  unburned 
gas  mixture  was  47  kPa  (~l/2  atm). 

Configuration  17.4.  Here  L  =  1187.8  cm,  d  =  17.4  cm,  which  is  similar  to  the  circular  cross-section  tube 
(diameter  17.4  cm)  used  in  [72],  Here,  the  right  end  of  the  channel  is  open  to  the  atmosphere,  and  the 
initial  gas  pressure  is  atmospheric,  as  in  the  experiments. 

Configuration  52.  Here  L  =  2130  cm,  d  =  52  cm,  which  is  similar  to  the  52  cm  diameter  tube  used  in 
[72],  where  the  channel  is  also  open  to  the  atmosphere  at  x  =  L. 


Table  3.  Configurations  for  Simulations 


Configuration 

7.6 

17.4 

52 

L  (cm) 

216.2 

1187.8 

2130 

d 

7.6 

17.4 

52 

£ 

1/3,  2/3 

0.3,  0.6 

0.3,  0.6 

x  =  l  boundary 

closed 

open 

open 

P0  (atm) 

0.464 

1 

1 

Cross-section,  experimental  channel 

square 

circular 

circular 

Experiment  reference 

[82] 

[72] 

[72] 

To  ignite  the  mixture  in  the  simulation,  we  placed  a  quarter-circular  region  of  hot,  burned  material  at 
the  left  wall  on  the  centerline  and  added  a  small  amount  of  extra  energy  to  the  burned  region.  The  additional 
energy  per  unit  mass  is  on  the  same  order  of  the  chemical  energy  release,  which  could  model  ignition  by  a 
low-energy  (~  lOOmJ)  spark.  The  resulting  shock  wave  is  very  weak  and  not  nearly  strong  enough  to  ignite 
a  detonation  directly.  It  does,  however,  cause  multiple  shock  reflections  and  shock-flame  interactions  that 
distort  and  wrinkle  the  flame  front. 

We  have  shown  that  it  is  possible  to  compute  the  behavior  of  a  propagating  laminar  flame  with  a 
resolution  of  5  computational  cells  across  the  thickness  of  the  reaction  front.  For  a  propagating  detonation, 
however,  we  required  over  200  computations  cells  across  the  detonation  front  thickness  to  properly  compute 
the  dynamics  of  the  unstable  detonation.  Here  we  are  most  concerned  with  the  acceleration  of  flame  and 
the  transition  to  detonation  and  less  interested  in  computing  the  large-scale  instabilities  of  the  detonation 
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front.  Hence,  in  the  simulations  described  below,  we  use  the  grid  resolution  required  for  the  laminar  flame 
thickness  and  do  not  try  to  fully  resolve  the  detonation  structure. 

3.1.3a.  Configuration  7.6 

In  the  experiments  [82],  obstacles  of  heights  1.27,  1.9,  and  2.53  cm  corresponded  to  blockage  ratios  of 
£  =  1/3,  1/2,  and  2/3,  respectively,  on  both  the  top  and  bottom  walls.  Both  ends  of  the  chamber  were 
closed.  The  pressure  of  the  stoichiometric  methane-air  mixture  inside  the  chamber  was  initially  47  kPa.  The 
mixture  was  ignited  by  an  electric  spark  at  the  centerline  of  the  channel. 

We  simulated  a  two-dimensional  rectangular  channel  (Figure  3.7)  with  obstacle  spacings  and  heights  as 
in  the  experiments.  After  ignition,  we  tracked  the  position  and  velocity  of  the  leading  edge  of  the  flame  front 
as  well  as  the  total  length  of  flame  surface  created  as  the  flame  evolves.  Flame  velocities  were  computed  at 
discrete  locations  along  the  length  of  the  channel  and  represent  an  average  velocity  over  the  interval  between 
two  successive  measurement  locations.  The  flame  surface  was  calculated  by  summing  the  total  length  of  the 
isosurface  on  which  Y  =  0.5  at  a  particular  instance  in  time.  Figure  3.8  compares  measured  and  computed 
flame  velocities  and  flame  surface  areas  for  /  =  1/3  and  £  =  2/3.  At  early  times  in  the  flame  development, 
the  simulations  and  experiments  showed  similar  flame  velocities  for  both  blockage  ratios.  Differences  arose 
further  downstream  (x  ~  100-150  cm)  as  the  flame  evolved. 

The  flame  acceleration  process  occurs  in  three  phases,  each  of  which  can  be  characterized  by  the  domi¬ 
nant  mechanism  driving  the  growth  in  flame  surface  area.  In  the  first  phase,  the  flame  is  folded  and  stretched 
by  a  laminar  flow  field  that  is  induced  by  the  thermal  expansion  of  the  combustion  products.  The  develop¬ 
ment  of  flame  surface  is  similar  for  /  =  1/3  and  2/3  for  this  phase,  x  <  70  cm  (Figure  3.8c).  Temperature 
maps  of  the  leading  edge  of  the  flame  as  it  passes  over  the  first  obstacle  are  shown  in  Figures  3.9a,b.  There 
is  relatively  little  wrinkling  of  the  flame  front  during  this  phase. 

In  the  second  phase,  the  predominant  mechanism  for  increasing  the  total  length  of  flame  surface  is  wrin¬ 
kling  by  fluid  dynamic  instabilities  (e.g.,  Kelvin-Helmholtz  and  Rayleigh-Taylor)  and  turbulent  fluctuations. 
Localized  regions  of  vorticity  stretch  and  fragment  the  flame  front,  thereby  increasing  the  total  amount  of 
flame  surface  area.  The  energy  released  at  the  flame  surface  causes  the  thermal  expansion  of  the  product 
gases,  which  causes  a  net  flow  through  the  channel.  Shear  layers  develop  downstream  of  obstacles  as  the 
fluid  is  accelerated  through  the  restricted  cross-sectional  area  above  them.  Fluid  dynamic  instabilities  and 
turbulence  in  the  shear  layers  contribute  to  generation  of  more  flame  surface  area.  The  flame  surface  and 
velocity  become  substantially  larger  for  £  =  2/3  than  £  =  1/3  for  x  >  70  cm.  In  this  case,  the  shear  layers 
develop  more  quickly,  since  the  flow  is  accelerated  to  a  higher  velocity  in  the  smaller  gap  between  obstacles. 

Temperature  maps  of  the  flame  fronts  near  x  =  100  cm,  Figures  3.9c  and  d,  show  that  there  is  much 
more  small-scale  flame  structure  and  therefore  much  more  flame  surface  in  the  £  =  2/3  case.  The  increase 
in  total  flame  surface  (see  Figure  3.8c)  continues  as  long  as  substantial  amounts  of  fuel  remain  in  between 
obstacles  behind  the  foremost  part  of  the  flame  front.  Then,  when  most  of  this  fuel  is  depleted,  the  flame 
surface  decreases  rapidly.  The  maximum  flame  surface  developed  for  £  =  2/3  is  larger  than  that  for  £  =  1/3 
because  the  increased  velocity  of  the  flame  front  allows  the  flame  to  propagate  farther  in  the  channel  before 
extended  reaction  zones  can  burn  out.  The  extra  amount  of  flame  surface  present  in  these  regions  results  in 
faster  depletion  of  the  fuel,  and  hence  the  steep  drop  in  flame  surface  for  x  >  130  cm. 

This  rapid  decline  in  flame  surface  is  slowed  as  the  flame  enters  the  third  phase.  When  the  speed  of  the 
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induced  flow  approaches  the  speed  of  sound  in  the  unburned  mixture  ahead  of  the  flame,  energy  released  at 
the  leading  edge  of  the  flame  front  generates  weak  pressure  waves  that  propagate  ahead  of  the  flame.  These 
pressure  waves  later  become  shocks,  which  can  reflect  from  obstacles  and  walls.  The  reflected  shocks  also 
collide  and  interact  with  portions  of  the  reaction  front.  The  flame  surface  is  wrinkled  by  these  shock-flame 
interactions  that  promote  Richtmyer-Meshkov  instabilities.  The  turbulence  generated  by  these  instabilities 
is  not,  necessarily,  homogeneous,  isotropic,  Kolmogorov  turbulence.  The  nature  of  this  turbulence  and  its 
interaction  with  the  flame  are  interesting  areas  of  future  investigation.  Here,  the  additional  energy  release 
caused  by  flame-surface  wrinkling  helps  to  sustain  the  flame  speed  and  slow  the  decline  in  the  net  flame 
surface.  This  process  occurs  for  x  >  160  cm  for  £  =  2/3,  but  does  not  begin  to  occur  in  the  £  =  1/3  case 
until  the  flame  reaches  the  end  of  the  domain.  A  longer  channel  would  be  necessary  to  observe  significant 
shock-flame  interactions  for  the  £  =  1/3  case.  Figures  3.9e,f  compare  the  temperature  maps  for  £  =  1/3  and 
2/3  when  the  leading  edge  of  the  flame  is  near  x  =  200  cm.  A  well-defined  shock  wave  has  formed  in  the 
£  =  2/3  case,  but  the  waves  ahead  of  the  flame  front  have  not  yet  coalesced  into  a  shock  for  £  =  1/3. 

3.1.3b.  Configurations  17.4  and  52 

In  the  larger-diameter  tubes  used  in  the  experiments  [72],  the  obstacles  were  annular  orifice  plates 
spaced  one  diameter  apart.  The  blockage  ratio  defined  in  these  experiments  is  £*  =  1  —  (D* /D')2,  where  D' 
and  D*  are  the  tube  diameter  and  orifice  diameter,  respectively.  One  end  of  the  tube  was  left  open  to  the 
atmosphere,  and  the  initial  gas  pressure  was  atmospheric  throughout  the  tube.  The  uniform  stoichiometric 
methane-air  mixture  was  ignited  near  the  tube  axis  at  the  closed  end.  Photodiodes  were  placed  at  various 
positions  along  the  walls  of  the  tubes,  and  reaction-front  velocities  were  calculated  based  on  time-of-arrival 
measurements . 

After  the  initial  acceleration  period,  two  regimes  of  propagation  velocity  were  found.  One  of  these, 
the  “choking”  regime  described  above,  is  characterized  by  a  velocity  close  to  1/2 Dcj-  The  second  regime, 
the  quasi-detonation  regime,  is  characterized  by  a  flame-front  velocity  just  less  than  Dqj ■  For  a  blockage 
ratio  of  l;*  =  0.3,  the  experimentally  measured  flame  speed  fluctuates  between  the  speeds  typical  of  these 
two  propagation  regimes.  This  indicates  that  the  17.4  cm  diameter  tube  is  close  to  the  critical  size  for 
detonation  propagation,  as  supported  by  the  observation  that  D'  / A  <  1  for  the  stoichiometric  methane-air 
system.  For  the  larger  blockage  ratio,  £*  =  0.6,  the  experimental  flame  velocity  approaches  a  steady  value 
of  approximately  700  m/s,  a  velocity  characteristic  of  the  choking  propagation  regime.  Similar  results  were 
obtained  for  the  52.0  cm  tube. 

We  performed  simulations  similar  to  these  experiments  using  the  configuration  shown  in  Figure  3.7  with 
d  =  D'  and  £  =  £*  for  configurations  17.4  cm  and  52  cm.  Because  of  the  differences  in  geometry,  the  obstacle 
heights  in  the  simulations,  h,  are  slightly  larger  than  the  heights  of  the  orifice  plates  in  the  experiments, 
h*  =  {D  —  D*)/ 2,  for  the  same  blockage  ratios.  Figure  3.10  shows  the  calculated  flame  velocities  and  surface 
areas  as  a  function  of  the  position  of  the  leading  edge  of  the  reaction  wave  for  £  =  0.3  and  0.6.  Figure 
3.11  shows  maps  of  temperature  (left  column),  pressure  (center  column),  and  velocity  (right  column)  for 
a  sequence  of  instances  over  the  course  of  the  simulation.  In  both  the  calculations  and  experiments,  the 
flame  accelerates  to  a  fixed  velocity  characteristic  of  the  choking  regime,  which  then  either  undergoes  DDT 
or  continues  to  propagate  at  this  average  speed.  The  initial  flame  acceleration  is  similar  to  that  described 
in  the  previous  section:  flames  are  stretched  by  the  thermal-expansion-induced  flow,  wrinkled  and  torn  by 
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turbulence  and  fluid  dynamic  instabilities,  and  further  fragmented  by  shock-flame  interactions.  For  example, 
the  black  line  in  each  frame  in  Figure  3.12  shows  the  progress  of  a  shock  colliding  with  and  passing  through 
a  flame.  As  the  shock  passes  through,  significantly  more  flame  surface  is  created  behind  it.  For  the  cases 
shown  in  Figure  3.10  and  3.11,  the  channel  is  long  enough  for  the  flames  to  progress  through  all  stages  of 
the  acceleration  process. 

The  evolution  of  flame  surfaces  for  0  <  x  <  450  cm  shown  in  Figure  3.10c  follows  the  same  trend  as 
in  the  early  stages  of  configuration  7.6.  For  x  >  450  cm  and  £  =  0.3,  the  flame  surface  sharply  decreases. 
At  this  point,  however,  the  reaction-front  velocity  jumps  to  Dqj,  indicating  that  a  detonation  was  initiated 
and  survived.  The  sequence  of  events  that  lead  to  DDT  is  shown  in  Figures  3.11i-m.  Strong  shock  waves 
formed  ahead  of  the  flame  front  reflect  off  of  the  channel  wall  and  the  faces  of  obstacles,  which  results  in  even 
stronger  waves  and  more  shock-flame  interactions.  Eventually,  Mach  stems  form  and  these  raise  the  local 
temperature  close  to  the  ignition  point.  These  regions  of  elevated  temperature,  or  hot  spots,  may  or  may  not 
ignite,  depending  on  the  ignition  delay  time  of  the  mixture  and  the  length  of  time  the  temperature  remains 
elevated.  The  hot  spot  created  by  a  Mach  stem  just  beginning  to  reflect  from  the  base  of  the  obstacle  in 
Figure  3. 1  lj  (rightmost  obstacle  in  figure)  leads  to  a  successful  detonation  ignition.  The  detonation  then 
propagates  into  unburned  fuel,  catches  up  to  the  leading  shock  wave,  and  proceeds  to  consume  nearly  all 
unburnt  fuel  in  that  region.  This  is  qualitatively  similar  to  the  process  reported  for  detonation  ignition  in 
hydrogen-air  mixtures  in  channels  with  obstacles  [6,  7].  The  size  of  the  system,  however,  is  considerably 
larger  for  the  methane-air  mixture. 

Velocity  curves  shown  (Figure  3.10a)  for  d  =  17.4  cm  and  £  =  0.3  indicate  that  the  detonation  propagates 
at  a  speed  much  less  than  Dqj-  This  reduction  in  velocity  is  a  result  of  recurring  detonation  diffractions  that 
continually  decouple  the  flame  from  the  leading  shock.  The  detonation  is  then  reignited  at  a  subsequent 
obstacle  in  the  manner  discussed  above.  Although  a  detonation  propagates  at  a  speed  greater  than  or 
equal  to  Dqj,  a  fast  deflagration  (decoupled  flame  and  shock)  propagates  significantly  slower.  Thus,  the 
time-averaged  velocity  for  this  quasi-detonation,  characterized  by  intermittent  periods  of  detonation  and 
fast  deflagration  propagation,  is  below  Dqj ■  An  example  of  detonation  failure  and  subsequent  reignition 
is  shown  in  Figures  3.11o-t.  This  repeated  ignition  and  decoupling  process  leads  to  the  observed  smaller 
propagation  velocities  for  this  case. 

Simulations  of  the  d  =  52  cm  channel  were  then  performed  for  £  =  0.3  and  £  =  0.6.  A  sequence  of 
temperature  (left  column),  pressure  (center  column),  and  velocity  (right  column)  maps  near  the  leading  edge 
of  the  flame  front  for  £  =  0.3  are  shown  in  Figure  3.13.  The  same  progression  of  events  observed  for  the 
d  =  17.4  cm  channel  can  be  seen.  The  first  occurrence  of  DDT  is  shown  in  panels  k  n.  Detonation  diffraction 
over  an  obstacle,  shock-flame  decoupling,  and  subsequent  DDT  is  shown  in  panels  o-t.  The  flame  velocities 
obtained  for  both  the  £  =  0.3  and  the  0.6  cases  are  shown  in  Figure  3.14.  For  £  =  0.3,  the  computed 
and  measured  flame  velocities  are  very  similar  during  the  initial  flame  acceleration  period  ( x  <  700  cm). 
Near  700  cm,  DDT  occurs  in  the  simulated  system  and  the  computed  velocity  jumps  to  ~  1800  m/s.  This 
propagation  speed  is  much  closer  to  Dqj  than  that  observed  in  configuration  17.4,  £  =  0.3  since  fewer 
instances  of  shock-flame  decoupling  take  place  in  the  larger  channel.  In  the  experiments,  DDT  first  occurs 
farther  downstream,  near  x  ~  1000  cm,  and  the  quasi-detonation  velocity  is  somewhat  smaller  than  the 
calculated  value.  For  £  =  0.6,  the  computed  flame  acceleration  is  close  to  the  experimental  data.  In  the 
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simulations,  several  instances  of  DDT  were  observed,  while  no  DDT  occurred  in  the  experiments.  The  total 
flame  surface  for  the  £  =  0.6  simulation  is  everywhere  greater  than  that  of  the  £  =  0.3  simulation.  Larger 
pockets  of  unburned  fuel  between  obstacles  take  longer  to  burn  and  delay  the  onset  of  the  rapid  decline  in 
flame  surface  area  that  occurred  in  configurations  7.6  and  17.4  for  the  larger  blockage  ratios.  By  the  time 
the  pockets  of  fuel  begin  to  burn  out,  the  leading  edge  of  the  propagating  flame  has  already  accelerated  to 
the  point  where  frequent  shock-flame  interactions  significantly  increase  the  amount  of  flame  surface,  leading 
to  less  rapid  net  losses  of  flame  surface  area. 

3.1.4.  Sensitivity  of  the  Model  to  Small  Parameter  Variations 

It  is  important  to  evaluate  the  sensitivity  of  the  computed  flame  acceleration  and  of  DDT  to  variations 
in  the  parameters  of  the  chemistry  model.  We  first  consider  the  sensitivity  of  the  flame  acceleration  by 
evaluating  how  much  velocity  profiles  change  when  there  are  systematic  variations  in  the  length  and  time 
scales  of  the  laminar  flame.  This  is  done  by  changing  model  parameters  A ,  q ,  Ea ,  and  7  to  create  moderate 
(10-15%)  variations  in  the  laminar  flame  velocity,  adiabatic  flame  temperature,  and  specific-heat  ratio.  We 
next  consider  the  impact  on  DDT  of  changing  these  parameters.  As  discussed  earlier,  model  parameters  also 
affect  various  length  and  time  scales  of  detonation  waves,  in  particular,  Dcj,  A,  and  ignition  delay  times. 


Table  4.  Configurations  for  Parameter  Sensitivity  Study 


V 

1 

2 

3 

4 

5 

A  (cm3g-s) 

1.64  x  1013 

1.2  x  1013 

1.64  x  1013 

4.411  x  1013 

1.64  x  1013 

Ea/  RTo 

67.55 

67.55 

69.45 

67.55 

67.55 

qM/RT0 

39.0 

39.0 

39.0 

34.71 

34.82 

7 

1.197 

1.197 

1.197 

1.197 

1.226 

Si  (cm/s) 

38.02 

32.0 

32.0 

38.02 

38.02 

Tb  (K) 

2210 

2210 

2210 

2000 

2210 

DCj  (nr/s) 

1820 

1820 

1820 

1724 

1854 

xd  (cm) 

0.229 

0.331 

0.356 

0.332 

0.0616 

The  first  question  to  address  is  how  changes  in  model  parameters  affect  the  velocity  of  the  flame  in  an 
obstructed  channel.  To  test  this,  we  performed  a  series  of  simulations  using  configuration  17.4  with  £  =  0.3 
and  0.6  and  varying  parameters  in  the  model.  A  description  of  the  parameter  sets  (denoted  V)  and  the 
corresponding  laminar  flame  and  ID  detonation  properties  is  given  in  Table  4. 

First  consider  different  sets  of  parameters  that  yield  the  same  T&  =  2210  K  but  different  values  of  the 
ID  laminar  flame  speed.  (In  Table  4,  these  are  VI-  V3.)  Figure  3.15a  shows  that  the  results  were  nearly 
identical  for  all  three  cases  with  £  =  0.6,  and  only  slight  differences  between  V3  and  VI  were  found  later 
in  the  flame  development  for  £  =  0.3.  These  simulations  suggest  that,  in  the  range  of  Si  considered,  the 
laminar  flame  speed  has  negligible  influence  on  the  evolution  of  the  flame  propagation  speed  in  the  channel. 

The  effect  of  the  adiabatic  flame  temperature  on  the  flame  acceleration  was  tested  by  changing  the 
amount  of  heat  released  per  unit  mass  of  fuel,  while  maintaining  a  constant  Si.  To  test  this,  we  compare  VI 
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and  W  Figure  3.15b  shows  the  flame  acceleration  computed  using  these  two  parameter  sets,  which  suggest 
this  200  K  temperature  difference  has  no  noticeable  effect  on  flame  acceleration.  The  effect  of  varying  7 
(and,  hence,  the  sound  speed)  on  the  flame  acceleration  is  shown  in  Figure  3.15c,  which  compares  VI  and 
V5.  Again,  only  slight  differences  between  the  two  cases  are  noticeable  for  both  £  =  0.3  and  0.6. 

The  second  important  issue  is  how  the  choice  of  model  parameters  affects  the  onset  of  the  DDT.  As 
we  discuss  later  in  this  report,  it  is  difficult  to  treat  this  issue  in  a  truly  quantitative  manner  because  DDT 
is  a  stochastic  process  that  depends  on  the  formation  of  relatively  small  hot  spots,  and  these  result  from 
combinations  of  shock  reflections  or  turbulent  fluctuations.  In  some  situations,  even  seemingly  imperceptible 
changes  in  any  physical  or  numerical  parameters  or  background  conditions  can  lead  to  significant  random 
variations  of  distances  or  times  to  detonation  initiation  [6] .  Even  with  these  caveats,  we  should  still  be  able 
to  get  a  qualitative  idea  of  the  likelihood  of  DDT  for  a  given  parameter  set. 

A  detailed  discussion  of  the  results  shown  in  Figure  3.16,  which  is  the  downstream  section  of  the  channels 
for  the  cases  shown  in  Figure  3.15,  is  given  in  [9].  Here  we  summarize  that  discussion  by  observing  that  the 
results  for  £  =  0.3  are  consistent  with  changes  in  chemical  induction  times  and  Ea,  which  in  turn  affect  the 
time  of  detonation  onset.  This  sensitivity  is  not  observed  for  £  =  0.6  which  does  not  show  the  transition. 
The  choice  of  7  (more  precisely,  7  —  1)  plays  a  large  role  in  detonation  initiation  and  propagation.  Using 
V5,  we  investigate  how  increasing  7  —  1  by  15%  (so  that  7  =  1.226)  affects  the  transition  to  detonation.  For 
£  =  0.3,  the  detonation  propagates  as  a  quasi-detonation  in  a  manner  similar  to  that  observed  using  set  VI. 
For  £  =  0.6  and  large  7,  there  are  several  instances  of  detonation  initiation  and  failure.  The  transitions  are 
infrequent,  and  the  average  propagation  velocity  is  smaller  than  Dcj  and  the  quasi-detonation  propagation 
velocities  observed  for  the  £  =  0.3  cases. 

The  results  shown  in  Figure  3.15  and  Figure  3.16  indicate  that  varying  the  model  parameters  to  result  in 
relatively  small  (10  —  15%)  changes  in  individual  laminar  flame  properties  has  little  impact  on  the  observed 
flame  acceleration.  In  general,  the  effects  on  DDT  are  also  small.  It  is  not  clear  whether  the  differences  in 
the  first  occurrence  of  DDT  among  the  several  parameter  sets  shown  in  Figure  3.16a  are  due  to  physical 
differences  in  the  modeled  systems  or  to  chance  fluctuations  in  the  thermodynamic  conditions  within  the 
hot  spots  that  initiate  detonations.  The  changes  in  model  parameters  for  V2  and  V3  both  have  the  effect 
of  increasing  induction  delay  times  behind  shocks  (an  effect  that  could  delay  DDT),  but  the  resulting  first 
occurrences  of  DDT  appeared  later  and  sooner  than  the  baseline  case,  respectively.  Small  changes  in  the 
thickness  of  the  detonation  wave  (and  detonation  cell  size)  result  in  correspondingly  small  changes  in  the 
average  propagation  velocity  and  first  occurrences  of  DDT  unless  the  system  is  near  a  critical  value  for 
detonation  propagation.  In  such  a  system,  even  a  small  increase  in  detonation  cell  size  could  impede  DDT. 
In  systems  sufficiently  larger  or  smaller  than  the  critical  size  for  DDT  for  a  given  fuel  mixture,  the  interaction 
of  the  induced  flow  and  shocks  with  the  wrinkled  flame  surface  seems  to  have  a  much  larger  influence  on  the 
large-scale  behavior  of  the  reaction  front  than  the  details  of  the  model  chemistry  in  these  types  of  obstructed 
channels. 

There  are  other  physical  properties  of  the  gas  mixture  that  have  not  been  discussed  here  that  could 
potentially  impact  the  first  occurrence  of  DDT.  For  instance,  the  ratio  of  the  acoustic  time  in  the  Zeldovich 
reactivity  gradient  to  the  reaction  time  [83]  has  been  shown  to  be  relevant  to  the  development  of  a  detonation 
[84].  It  is  not  possible  to  vary  the  reaction  time  and  the  induction  time  independently  using  a  one-step 
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reaction  model.  This  could  be  done  with  a  multiple-step  reaction  model,  but  at  a  higher  computational  cost. 
Evaluating  whether  solutions  obtained  using  a  multiple-step  reaction  model  would  produce  more  accurate 
results  and,  if  so,  whether  the  improvement  justifies  the  additional  computational  costs  is  left  for  future 
investigations.  Likewise,  we  have  not  discussed  how  changes  in  the  laminar  flame  thickness  affect  flame 
acceleration.  The  thickness  of  the  flame  influences  its  ability  to  become  wrinkled  and  distorted  by  shocks 
and  fluid-dynamic  instabilities,  which  in  turn  greatly  affects  the  total  energy  release  rate.  We  will  return  to 
this  important  issue  in  Sections  3.3  and  3.4. 
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(a)  (b) 


Figure  3.1.  Parametric  curves  for  which  (a)  Tj,  =  2210  K,  Dqj  =  1820  m/s  and  (b)  Si  =  38.02  cm/s, 
Xd  =  0.229  cm.  The  points  of  intersection  in  the  two  figures  give  the  values  of  qM/RT0 ,  7,  Ea/RT0 ,  and  A 
used  in  conjunction  with  the  single-component  reaction  model  (Eq.  3.3). 


x  (cm) 


Figure  3.2.  Temperature  (red)  and  reaction-rate  (green)  profiles  for  a  laminar  flame  calculated  using 
Eqs.  (3. 5-3. 6)  (lines),  two-dimensional  Navier-Stokes  equations  (Eqns.  2. 1-2. 5)  with  dxmin  =  0.018125  cm 
(triangles),  and  the  two-dimensional  Navier-Stokes  equations  (Eqns.  2. 1-2. 5)  with  dxmin  =  0.00113  cm 
(squares) . 
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Figure  3.3.  Front-propagation  velocity  as  a  function  of  time  after  the  denotation  has  reached  a  quasi-steady 
propagation  speed  of  1813  m/s  (solid  line)  and  Dcj  for  a  stoichiometric  (9.5%)  nrethane-air  mixture  (dashed 
line). 
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Figure  3.4.  Selected  density  maps  of  the  full  channel  width  in  the  vicinity  of  the  detonation  front  at  times 
(a)  t  =  2.480  ms,  (b)  t  =  2.505  ms,  (c)  t  =  2.542  ms,  (d)  t  =  2.586  ms,  (e)  t  =  2.620  ms,  (f)  t  =  2.648  ms, 
and  (g)  t  =  2.677  ms  after  initiation. 
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Figure  3.5.  Temperature  maps  near  (a)  a  triple-shock  configuration  to  which  a  transverse  detonation 
is  attached,  (b)  a  triple-shock  configuration  where  detonation  occurs  at  the  leading  edge  but  not  in  the 
transverse  wave,  and  (c)  a  triple-shock  structure  showing  characteristics  of  both  (a)  and  (b),  and  for  which  the 
transverse  wave  is  a  detonation  between  the  oblique  shock  and  slip  line  and  an  inert  shock  farther  away  from 
the  leading  edge.  Structures  shown  in  the  figure  include  (1)  leading  inert  shocks,  (2)  transverse  detonations, 
(3)  induction  zones  (unreacted  material),  (4)  oblique  shocks,  (5)  slip  lines,  (6)  overdriven  detonations,  (7) 
decoupled  shocks  and  reaction  zones,  (8)  pockets  of  unreacted  material,  and  (9)  inert  transverse  shocks. 


45 


(a) 


t  =  2.JV0  ms 


(C) 


x  (cm) 

I  =  2.645  ms 


x  (cm) 


x  (cm) 


I  =  2.66  n 


- C - 

ym 

M 

V 

_ 

\ 

] 

t  =  2.547  ms 


i  =  2.554  ms 


x  (cm) 


Density  (g/cm5) 


0.0005  O.COtS  00055  0.0055  00045  0.0055  0.0065  0.007 


x  (ctn) 


x(cm) 


Figure  3.6.  Density  maps  near  several  representative  triple  point  collisions:  (a)  two  triple  points  with  inert 
transverse  waves,  (b)  one  triple  point  with  transverse  detonation  and  one  with  an  inert  transverse  wave,  and 
(c)  two  triple  points  with  transverse  detonations.  Time  increases  from  left  to  right  in  each  figure. 
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Figure  3.7.  Schematic  of  general  two-dimensional  channel  with  obstacles  used  in  DDT  calculations.  We  as¬ 
sume  the  channel  is  symmetric  and  simulate  only  the  lower  half,  so  symmetry  conditions  ( du/dy  =  dY/dy  = 
dT/dy  =  dP/dy  =  v  =  0)  are  applied  at  the  channel  center  line,  y  =  d/2.  The  obstacles  are  taken  to  be  2  cm 
thick,  and  their  heights  are  set  based  on  the  desired  blockage  ratio,  £  =  2 h/d.  Spacings  between  obstacles 
are  set  equal  to  d,  and  the  obstacles  fill  the  entire  length  of  the  channel.  Walls  and  obstacle  surfaces  are 
adiabatic  no-slip  reflecting  boundaries.  Calculations  are  initialized  by  placing  a  quarter-circular  region  of 
burned  material  in  the  upper  left  corner  of  the  domain.  This  initial  flame  radius  is  0.25  cm  in  all  calculations. 


(a)  (b) 

§=1/3  §=2/3 


X  (m) 

Figure  3.8.  Configuration  7.6:  (a  and  b)  Computed  and  measured  [82]  flame-propagation  velocities  and 
(c)  computed  flame-surface  as  a  function  of  the  position  of  the  leading  edge  of  the  reaction  front. 
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Figure  3.9.  Temperature  maps  near  the  leading  edge  of  the  flame  for  configuration  7.6  with  £  =  1/3  (left) 
and  £  =  2/3  (right)  at  several  locations  throughout  the  channel:  (a  and  b)  near  first  obstacle,  (c  and  d)  near 
x  =  100  cm,  and  (e  and  f)  near  x  =  200  cm.  The  top  temperature  scale  is  for  burned  material,  and  the 
bottom  scale  is  for  unburned  material.  Time  increases  from  top  to  bottom  in  each  column. 


(a)  (b) 


Figure  3.10.  Configuration  17.4:  (a  and  b)  Computed  and  measured  [72]  reaction  front-propagation  veloc¬ 
ities  and  (c)  calculated  flame-surface  lengths  as  a  function  of  the  position  of  the  leading  edge  of  the  reaction 
front. 
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Figure  3.11.  Legends  and  caption  on  next  page. 
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Figure  3.11.  Configuration  17.4:  Temperature,  pressure,  and  x-velocity  sequences  (see  previous  page)  in 
the  simulation  of  flame  acceleration  and  DDT  in  an  obstructed  channel  with  £  =  0.3.  Pressure  and  velocity 
scales  change  with  time  as  the  static  pressure  and  bulk  velocity  in  the  channel  grow  while  the  flame  is 
accelerating.  There  are  two  temperature  scales,  as  indicated  in  the  legends  above.  The  upper  scale  (cooler) 
refers  to  the  unburned  gas,  the  lower  scale  (hotter)  refers  to  the  burned  gas.  The  letters  on  the  left  side  refer 
to  the  frames  on  the  previous  page. 
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Figure  3.12.  Temperature  maps  near  the  leading  edge  of  flame  in  configuration  17.4  with  £  =  0.3  as  a 
shock  interacts  with  the  flame  front.  The  heavy  black  line  indicates  the  location  of  a  shock,  and  the  arrow 
indicates  the  direction  of  propagation.  The  top  temperature  scale  is  for  burned  material,  and  the  bottom 
scale  is  for  unburned  material.  Time  increases  from  top  to  bottom. 
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Figure  3.13.  Legends  and  caption  on  next  page. 
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Figure  3.13.  Configuration  52:  Temperature,  pressure,  and  x- velocity  sequences  (see  previous  page)  in  the 
simulation  of  flame  acceleration  and  DDT  in  an  obstructed  channel  with  £  =  0.3.  Pressure  and  velocity  scales 
change  with  time  as  the  static  pressure  and  bulk  velocity  in  the  channel  grow  while  the  flame  is  accelerating. 
There  are  two  temperature  scales,  as  indicated  in  the  legends  above.  The  upper  scale  (cooler)  refers  to  the 
unburned  gas,  the  lower  scale  (hotter)  refers  to  the  burned  gas.  The  letters  on  the  left  side  refer  to  the 
frames  on  the  previous  page. 
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(a)  (b) 
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Figure  3.14.  Configuration  52:  (a  and  b)  Computed  and  measured  [72]  reaction  front-propagation  velocities 
and  (c)  calculated  flame-surface  lengths  as  a  function  of  the  position  of  the  leading  edge  of  the  reaction  front. 
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Figure  3.15.  Reaction  front-propagation  velocities  as  a  function  of  the  position  of  the  leading  edge  of  the 
reaction  front  for  x  <  300  cm  in  configuration  17.4  calculated  using  parameter  sets,  V ,  (a)  1,  2,  3,  (b)  1,  4, 
and  (c)  1,  5.  See  Table  4  for  descriptions  of  the  parameter  sets.  Lines  with  symbols  represent  £  =  0.6,  and 
those  without  represent  £  =  0.3. 
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Figure  3.16.  Reaction  front-propagation  velocities  as  a  function  of  the  position  of  the  leading  edge  of  the 
reaction  front  for  x  >  300  cm  in  configuration  17.4  calculated  using  parameter  sets,  V ,  (a)  1,  2,  3,  (b)  1,  4, 
and  (c)  1,  5.  See  Table  4  for  descriptions  of  the  parameter  sets.  Lines  with  symbols  represent  £  =  0.6,  and 
those  without  represent  £  =  0.3. 
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3.2.  Stochasticity  and  Resolution  (Task  4) 

Stochasticity  is  the  ability  of  a  process  to  deviate  randomly  from  its  mean  path.  Experimental  observations 
show  that  DDT  is  a  stochastic  process  in  the  sense  that  there  is  some  uncertainty  in  time  and  location  for  the 
detonation  initiation.  This  uncertainty  is  an  important  aspect  of  the  DDT  that  involves  multiple  stochastic 
phenomena,  including  flow  instabilities,  turbulence,  many  interactions  between  shocks,  flames,  and  vortices, 
and  the  resulting  hot-spot  formation.  Since  the  detonation  usually  arises  from  one  of  many  hot  spots  that 
stochastically  appear  in  the  system,  detonation  initiation  also  becomes  stochastic.  Small  fluctuations  in 
density,  temperature,  and  composition  play  an  important  role  in  the  development  of  stochasticity  in  experi¬ 
mental  systems,  which  are  also  affected  by  uncertainties  in  initial  conditions.  The  real  cause  of  stochasticity, 
however,  is  embedded  in  the  complexity  of  underlying  physical  phenomena  and  is  present  in  RNSE-based 
numerical  models  that  describe  these  phenomena. 

High-Reynolds  number  flows,  those  for  which  we  can  only  reasonably  resolve  a  portion  of  the  range  of 
turbulent  scales,  are  stochastic  in  the  sense  that  they  are  very  sensitive  to  minute  changes  in  computational 
parameters.  Small  variations  in  reaction  rates,  initial  conditions,  background  conditions,  features  of  the 
computational  grid,  and  even  in  the  numerical  roundoff  that  occurs  at  a  particular  time  in  the  simulation, 
can  change  the  exact  time  and  location  when  various  features  of  the  flow  arise  and  interact  or  are  quenched. 
The  issue  of  stochasticity  in  computations  of  DDT  is  closely  related  to  numerical  resolution  tests,  since 
changing  the  computational  cell  size  affects  a  non-converged  solution,  and  triggers  a  stochastic  response  at 
the  same  time.  When  stochastic  deviations  are  comparable  to  the  systematic  effect  of  numerical  resolution, 
the  effects  of  resolution  cannot  be  properly  evaluated  without  estimating  the  dispersion  of  the  results  related 
to  stochastic  properties  of  the  system. 

How  such  small-scale  perturbations  affect  the  large-scale  behavior  of  a  system  is  a  question  largely 
unaddressed  in  turbulent  reacting  flows  because  of  their  physical  complexity  and  the  high  costs  of  the 
simulations.  It  is  important,  however,  to  understand  whether  a  particular  DDT  mechanism  is,  in  fact, 
robust  to  small  variations  in  a  turbulent,  stochastic  flow. 

3.2.1.  Background 

Our  previous  experience  with  DDT  simulations  has  shown  that  some  systems  are  more  susceptible  to  the 
stochastic  behavior  than  others,  depending  on  geometries  and  initial  conditions  involved. 

In  the  earliest  DDT  computations  for  low-pressure  acetylene  and  ethylene  [41,43,  44],  we  did  not  observe 
any  significant  stochastic  behavior,  though  no  systematic  studies  of  stochasticity  were  performed  for  these 
systems.  These  computations  showed  excellent  convergence  with  changes  in  the  computational  grid.  This 
suggests  that  systematic  effects  of  numerical  resolution  were  not  altered  by  any  significant  random  devia¬ 
tions.  Stochastic  effects  were  observed  in  computations  of  DDT  in  obstructed  channels  [6]  that  included 
both  numerical  convergence  studies  and  dedicated  stochasticity  tests.  The  latter  were  based  on  background 
temperature  variations  in  the  range  dT  =  ±0.01K  that  were  too  small  to  cause  any  systematic  effects.  All 
other  model  parameters,  including  the  numerical  resolution,  were  kept  constant.  Multiple  computations  per¬ 
formed  for  these  slightly  different  background  temperatures  showed  practically  no  effect  on  the  initial  flame 
acceleration,  but  revealed  moderate  stochastic  dispersion  for  the  run-up  distance  to  DDT.  This  dispersion 
was  of  the  order  of  20%  for  the  larger  obstacle  spacing  S  >  d  (Configuration  1  in  Section  2.3.3),  and  of 
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the  order  of  40%  for  the  smaller  obstacle  spacing  S  <  d  (Configuration  2  in  Section  2.3.3).  These  tests 
were  repeated  for  several  numerical  resolutions  and  showed  that  stochastic  deviations  were  comparable  to 
systematic  effects  of  numerical  resolution. 

3.2.2.  Methane  Tests 

Now  we  describe  the  stocliasticity,  model  sensitivity,  and  numerical  resolution  tests  performed  for  flame 
acceleration  and  DDT  in  the  stoichiometric  methane-air  mixture.  We  first  describe  these  tests  individually, 
and  then  discuss  the  implications  for  producing  a  predictive  computational  ability  for  DDT. 

3.2.2a.  Laminar  Flame  Resolution 

The  numerical  issue  of  how  much  grid  resolution  is  necessary  at  the  reaction  front  becomes  extremely 
important  when  attempting  to  simulate  flame  acceleration  and  DDT  for  very  large-scale  explosions.  A 
major  constraint  on  the  size  of  the  system  that  can  be  simulated  is  set  by  the  number  of  small  cells  required 
to  obtain  a  sufficiently  accurate  description  of  the  reacting  flow.  Thus,  it  is  useful  to  know  the  largest  cell 
size  that  can  be  used  to  meet  this  requirement.  Our  previous  tests  have  shown  that  the  amount  of  resolution 
required  to  produce  the  correct  laminar  flame  speed  depends  strongly  on  the  complexity  of  kinetics  and 
diffusion  used  in  the  model.  We  have  seen  above  that  the  methane-air  flame  with  a  single  Arrhenius  step, 
which  is  a  rather  simplified  chemical-diffusion  model,  needs  to  be  resolved  with  only  three  cells  to  propagate 
with  the  approximately  correct  velocity.  In  previous  studies  [85],  which  used  full  chemical  reaction  models, 
we  found  that  we  needed  20-50  cells  to  resolve  the  flame  thickness,  depending  on  the  number  of  species 
in  the  model.  Because  of  the  significant  sensitivity  of  the  flame  speed  and  subsequent  acceleration  for 
multispecies  models,  resolution  tests  of  such  models  were  even  more  crucial  if  the  results  were  expected  to 
be  quantitatively  or  even  qualitatively  correct. 

The  next  level  question,  one  that  goes  beyond  isolated  laminar-flame  resolution,  is  what  is  the  resolution 
required  to  describe  the  flame  correctly  in  the  flame-acceleration  region  of  an  entire  flow  calculation.  This 
depends  not  only  on  the  resolution  at  the  front,  but  the  influence  of  other  properties  of  the  flows,  including 
instabilities  and  turbulence  that  also  depend  on  resolution.  Here  we  break  this  issue  into  two  parts.  First, 
we  have  addressed  the  question:  What  is  the  resolution  needed  at  the  front  to  achieve  the  correct  flame 
acceleration? 

We  tested  the  effect  of  the  reaction  front  resolution  on  flame  acceleration  in  configuration  7.6  defined 
in  Table  3  with  £  =  1/3  using  the  model  parameters  shown  in  Table  2.  In  all  cases,  the  flame  was  initiated 
from  a  spark  of  size  rj  =  0.25  cm.  Three  values  of  dxmin  were  tested,  0.018124  cm  (corresponding  to  ~  3 
computational  cells  per  laminar  flame  thickness),  0.009062  cm  (~  5  cells/flame  thickness),  and  0.004531  cm 
(~  10  cells/flame  thickness).  Figure  3.17a  shows  the  computed  propagation  velocity  of  the  leading  edge  of 
the  flame,  and  Figure  3.17b  shows  the  development  of  the  flame  surface,  both  as  a  function  of  the  position 
of  the  leading  edge  of  the  flame.  As  the  minimum  grid  size  is  decreased,  and  the  flame  front  becomes  more 
resolved,  the  peak  amount  of  flame  surface  increases.  The  computed  flame  velocity  during  the  later  stages  of 
the  flame  acceleration  process  also  increases  when  the  number  of  grid  cells  in  the  reaction  zone  is  increased 
from  3  to  5.  A  subsequent  increase  from  5  cells  to  10  cells  has  little  impact  on  the  bulk  acceleration  of  the 
flame. 
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3.2.2b.  Background  Flow  Resolution 

Our  numerical  model  includes  an  adaptive  mesh  refinement,  which  means  that  the  computational  cell 
size  in  the  vicinity  of  flames  and  shocks  is  significantly  smaller  than  in  the  background  flow,  away  from  steep 
gradients.  Since  the  maximum  and  minimum  computational  cell  sizes  can  be  changed  independently,  we  can 
consider  the  effect  of  the  resolution  of  the  background  flow  on  flame  acceleration.  There  are  two  case  studies 
for  this:  one  for  configuration  7.6  and  another  for  configuration  17.4  (Table  3).  First,  consider  configuration 
7.6,  which  is  shown  in  Figure  3.18.  Here,  the  size  of  the  finest  computational  cells  in  the  vicinity  of  shocks 
and  flame  fronts  were  held  fixed  at  dxmin  =  0.018124  cm  while  the  largest  grid  sizes  in  the  background  flow 
were  systematically  decreased  from  dxmax  =  a  =  0.211  cm  to  a/8  =  0.0264  cm.  The  flame  acceleration 
increases  noticeably  as  the  grid  sizes  are  decreased  from  a  to  a/2  to  a/4.  Only  small  differences  between  the 
curves  for  a/4  and  a/8  are  visible,  suggesting  some  level  of  grid  convergence  has  been  attained. 

Next,  consider  configuration  17.4,  which  is  the  physically  more  difficult  case  because  it  leads  to  DDT, 
and  the  simulations  summarized  by  the  curves  in  Figure  3.19.  As  in  the  previous  figure,  the  three  curves 
represent  cases  for  which  we  apply  the  same  initial  physical  conditions  and  finest  grid  size,  but  vary  the 
size  of  the  coarsest  computational  grids.  The  two  more  resolved  cases  give  very  similar  flame  acceleration 
until  just  before  DDT.  In  the  least  resolved  case,  the  flame  acceleration  is  notably  slower.  We  can  therefore 
consider  the  two  more  resolved  cases  as  having  acceptable  background  resolution.  One  explanation  for  this 
variation  in  flame  acceleration  from  coarse  to  fine  grids  is  that  for  finer  background  resolution,  more  of 
the  turbulent  length  scales  are  resolved.  This  increase  in  resolution  of  the  turbulence,  in  turn,  means  that 
finer-scale  background  effects  contribute  to  flame-surface  wrinkling,  there  is  faster  energy  release,  and  hence 
faster  flame  acceleration. 

3.2.2c.  Stochasticity,  Resolution,  and  DDT 

The  issue  of  multiple,  interacting  stochastic  processes  in  DDT  is  one  of  the  most  difficult  and  important 
issues  for  safety  and  risk  analysis.  How  much  variation  in  the  distance  to  DDT  can  we  expect,  given  the 
stochastic  nature  of  turbulence  and  hot-spot  formation? 

To  investigate  this  point,  consider  again  configuration  17.4.  The  two  curves  representing  the  most 
resolved  simulations  from  Figure  3.19  are  reproduced  in  Figure  3.20.  The  results  of  four  other  nearly  identical 
simulations  are  also  shown.  The  three  green  curves  represent  simulations  performed  with  dxmax  =  0.145  cm, 
and  for  the  three  blue  curves,  dxrnax  =  0.0725  cm.  The  differences  between  the  curves  of  each  particular 
color  were  minute  variations  in  the  initial  system  temperatures  (0.01  K).  These  simulations,  which  are 
nearly  identical  physically,  show  a  range  of  locations  and  timings  in  which  DDT  occurred.  In  fact,  all  of 
the  simulations  shown  in  Figure  3.20  have  undergone  DDT  at  about  the  same  flame  and  shock  velocity, 
~1000  m/s.  In  all  cases,  the  flame  acceleration  is  similar,  but  the  transition  occurs  anywhere  in  about  a 
50  cm  region,  which,  equates  to  three  or  four  obstacle  spacings.  This  corresponds  to  uncertainty  of  about 
15%  in  the  run-up  distance  to  DDT,  which  is  consistent  with  20%  uncertainty  obtained  in  similar  tests  for 
lrydrogen-air  mixtures  discussed  above. 
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(a) 


(b) 


Figure  3.17.  (a)  Flame  velocity  and  (b)  flame-surface  length  as  a  function  of  position  of  leading  edge  of 
a  flame  for  configuration  7.6,  £  =  1/3  using  three  successively  higher  grid  resolutions  at  the  flame  front, 
representing  3,  5,  and  10  grid  cells  per  laminar  flame  thickness. 


x  (m) 

Figure  3.18.  Flame  velocity  as  a  function  of  position  of  leading  edge  of  a  flame  for  configuration  7.6,  £  =  1/3 
using  three  successively  higher  grid  resolutions  in  the  background  flow ,  so  that  the  largest  computational  cells 
in  the  calculations  were  successively  decreased  by  factors  of  2,  ranging  from  dxmax  =  0.211  cm  to  0.026376  cm. 
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Figure  3.19.  Flame  velocity  as  a  function  of  position  of  leading  edge  of  a  flame  for  configuration  17.4,  £  =  0.3 
using  three  successively  higher  grid  resolutions  in  the  background  flow ,  so  that  the  largest  computational  cells 
in  the  calculations  were  successively  decreased  by  factors  of  2,  ranging  from  dxmax  =  0.29  cm  to  0.07256  cm. 


x  (m) 


Figure  3.20.  Flame  velocity  as  a  function  of  position  of  leading  edge  of  a  flame  for  configuration  17.4, 
£  =  0.3  for  six  nearly  identical  systems.  Green  lines  represent  calculations  done  with  dxmax  =  0.145  cm. 
Blue  lines  represent  calculations  done  with  dxmax  =  0.07256  cm.  For  each  color,  different  lines  represent 
different  initial  system  temperatures:  297.99  K  (solid  lines),  298.00  K  (dashed  lines),  298.01  K  (dash-dotted 
lines).  Onset  of  DDT  is  found  to  vary  stochastically  for  these  numerical  and  background  perturbations. 
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3.3.  DDT  in  Fuel-Lean  and  Fuel-Rich  Mixtures  (Task  2) 

Above  we  outlined  a  systematic  approach  to  developing  a  chemical-diffusive  submodel  for  stoichiometric 
methane-air  mixtures.  Based  on  the  success  of  this  model  when  it  was  used  in  full  numerical  simulations  for 
predicting  flame  acceleration  and  DDT,  we  proceeded  using  a  similar  approach  to  develop  models  for  both 
lean  and  rich  methane-air  mixtures.  The  ultimate  objective  then  is  to  produce  some  model  that  could  be 
interpolated  when  used  for  gradients  of  stoichiometry,  which  is  what  we  normally  would  see  in  a  naturally 
developing  mixture  of  gases.  Below  we  describe  a  reasonably  good  model  for  lean  mixtures  of  methane  and 
air  and  show  that  this  approach  does  not  generalize  directly  to  rich  mixtures.  Several  alternate,  more  general 
approaches  for  the  form  of  a  chemical-diffusive  model  are  discussed.  The  key  condition  for  computing  DDT 
is  that  the  same  model  must  work  reasonably  well  for  both  flames  and  detonations.  Figure  3.21  shows  how 
several  important  properties  of  methane-air  mixtures  vary  with  stoichiometry  [33,  57,  58,  73]. 

3.3.1.  Fuel-Lean  Model 

In  order  to  find  a  model  that  works  well  in  full-scale  simulations,  we  look  for  overlap  in  the  sets  of  model 
parameters  that  will  give  reasonable  values  for  both  flames  and  detonation.  Calibration  is  performed  using 
the  same  methodology  described  in  Section  3.1.1.  The  target  values  for  a  fuel-lean,  7%  nrethane-air  system 
are  given  along  with  those  for  a  stoichiometric  mixture  in  Table  1.  Figure  3.22  shows  results  for  key  output 
parameters  as  a  function  of  the  input  variables.  As  before,  we  select  values  from  intersections  of  these  curves. 
For  this  lean  case  we  find  both  7  and  q  are  slightly  smaller  than  those  obtained  for  the  stoichiometric  case  as 
one  might  expect  since  the  total  amount  of  fuel  and,  hence,  energy  that  can  be  released  is  also  smaller.  There 
is  considerable  uncertainty  in  the  detonation  cell  size  for  this  mixture,  so  two  curves  are  shown  representing 
loose  bounds  on  the  acceptable  size  of  the  computed  Xd-  This  figure  gives  a  range  of  values  A,  Ea  that  lie 
along  the  Si  =  23.7  cm/s  curve  bounded  by  the  two  dashed  curves.  We  choose  a  value,  Ea  =  57.55 RT0,  that 
lies  in  the  middle  of  this  range.  The  corresponding  value  of  A  and  all  other  calibrated  input  parameters  are 
given  in  Table  5. 

We  next  consider  the  behavior  of  a  flame  calculated  using  the  set  of  parameters  in  Table  5  for  configu¬ 
ration  17.4.  Here  we  only  investigate  the  case  for  which  £  =  0.3.  The  resulting  flame  velocity  is  shown  as  a 
function  of  position  in  the  channel  in  Figure  3.23.  Flame  acceleration  is  nominally  slower  than  that  observed 
for  the  stoichiometric  case  discussed  in  Section  3.1.3,  consistent  with  the  fact  that  Si  is  smaller  and  Xi  is 
slightly  larger  for  the  lean  flame  than  they  are  for  stoichiometric  flames.  It  is  worth  noting  here  that  the 
computed  xi  =  0.067  cm  is  still  about  a  factor  of  three  smaller  than  measured  values  [33],  which  could  have 
an  impact  on  the  evolution  of  flame  surface  and,  hence,  acceleration.  Experimental  flame  velocity  measure¬ 
ments  [72]  for  6.5%  and  8%  mixtures  are  also  shown  in  Figure  3.23.  These  curves  bound  the  velocity  profile 
of  the  simulated  7%  mixture,  suggesting  that  the  model  parameters  given  in  Table  5  at  least  qualitatively 
predict  the  right  flame  acceleration.  Figure  3.23  also  shows  the  results  of  two  parametric  sensitivity  analyses 
for  which  T/  and  Si  were  systematically  varied.  The  results  are  the  same  as  those  discussed  in  Section  3.1.4, 
namely  that  small  changes  in  these  flame  properties  have  little  effect  on  the  overall  flame  acceleration.  Also, 
no  DDT  was  observed  for  this  mixture.  In  the  experiments,  an  8%  mixture  was  found  to  undergo  DDT  while 
a  6.5%  mixture  did  not. 
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Table  5.  Material,  chemical,  and  reaction-wave  parameters  for  lean  (7%)  methane-air 


Quantity 

Value 

Definition 

Input 

Po 

1  atm 

Initial  pressure 

T0 

298  K 

Initial  temperature 

Po 

1.1042  x  10~3  g/cm3 

Initial  density 

7 

1.18 

Adiabatic  index 

Ms 

27 

Molecular  weight 

A 

3.2  x  1012  cm3/g/s 

Pre-exponential  factor 

Ea 

57.55  R  To 

Activation  energy 

q 

36.8  R  To  /M 

Chemical  energy  release 

vo 

3.6  x  10~6  (g/s-cm-K0-7) 

Transport  constants 

7i 

o 

II 

to 

o 

6.25  x  10"6  (g/s-cm-K0-7) 

Transport  constants 

Output 

Si 

23.7  cm/s 

Laminar  flame  speed 

Tb 

1971  K 

Adiabatic  flame  temperature 

Xl 

0.067  cm 

Laminar  flame  thickness 

Dcj 

1692  m/s 

Chapman-Jouguet  detonation  velocity 

Xd 

1.01  mm 

ID  half-reaction  detonation  thickness 

3.3.2.  Fuel-Rich  Model 

Finding  an  adequate  model  of  the  form  we  are  using  for  fuel-rich  mixtures  is  much  more  difficult.  The 
fundamental  problem,  shown  in  Figure  3.21,  is  that  for  methane- rich  mixtures  the  rate  of  change  of  detonation 
properties  (e.g.  A  and  Dqj)  as  a  function  of  methane  concentration  is  qualitatively  different  from  that  of 
methane-lean  mixtures.  The  detonation  velocity  and  cell  size  remain  larger  and  smaller,  respectively,  for 
higher  methane  concentrations  than  do  the  laminar  flame  velocity  and  adiabatic  flame  temperature.  Because 
of  this,  curves  of  constant  S';  and  Xd  for  this  particular  mixture  do  not  intersect  within  a  range  of  realistic 
values  of  the  activation  energy  as  shown  in  Figure  3.24.  The  figure  shows  curves  for  two  different  S;  that 
reflect  the  large  variation  in  measured  flame  speeds  found  in  the  literature  and  two  curves  for  detonation 
thickness  to  reflect  the  uncertainty  in  actual  detonation  cell  size  and  the  correlation  factor  of  the  cell  size  to 
Xd-  Still,  there  is  no  intersection.  Here,  we  use  two  different  approaches  to  address  this  difficulty. 

In  the  first  approach,  we  strike  a  compromise  that  attempts  to  match  flame  parameters  more  closely 
than  detonation  properties.  The  justification  for  this  is  that  flame  acceleration  and  onset  of  DDT  depend 
more  strongly  on  the  properties  of  the  flame  than  on  the  properties  of  a  detonation  wave.  These  are  given 
in  the  first  column  of  Table  6.  The  resulting  adiabatic  flame  temperature  and  laminar  flame  speed  lie  just 
outside  the  ranges  given  in  Table  1,  while  the  detonation  velocity  and  the  detonation  cells  are  noticeably 
smaller  than  the  ranges  given  in  Table  1.  The  laminar  flame  thickness  obtained  using  these  parameters 
(. Xi  =  0.0631  cm),  while  larger  than  that  obtained  for  a  stoichiometric  mixture,  is  significantly  smaller  than 
that  measured  in  experiments  [33]. 
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In  the  second  approach,  we  attempt  to  bring  the  bounds  of  the  flame  speed  curves  within  the  bounds  of 
the  detonation  thickness  curves  by  choosing  a  smaller  value  of  A  for  each  particular  Ea  and  then  varying  the 
thermal  conductivity,  a  parameter  we  had  assumed  to  be  fixed  and  equal  to  that  of  air  in  our  stoichiometric 
and  lean  models,  to  speed  up  and  artificially  thicken  the  laminar  flame.  This  procedure  has  no  effect  on 
the  detonation  wave  since  thermal  conduction  does  not  play  a  significant  role  in  that  combustion  mode. 
Figure  3.25  shows  four  different  curves  in  the  n0-Ea  parameter  space.  These  curves  represent  the  values 
of  k0  needed  to  achieve  either  Si  =  15.7  cm/s  or  23.7  cm/s  for  the  values  of  A  shown  in  Figure  3.24  that 
define  the  bounds  Xd  =  0.69  cm  and  1.1  cm.  Hence,  the  top  and  bottom  curves  in  the  figure  define  the 
bounds  on  the  values  of  kq  that  can  be  used  to  obtain  Si  and  Xd  within  the  ranges  shown  in  Table  1.  We 
are  thus,  again,  in  the  situation  where  the  choice  of  parameters  is  non-unique.  The  additional  dimension  in 
the  parameter  space,  now  occupied  by  k,  allows  us  to  vary  the  laminar  flame  thickness  independently  of  the 
flame  speed  in  an  attempt  to  create  a  flame  that  is  more  physically  realistic.  For  fixed  A1  Ea,  decreasing  kq 
results  in  a  reduction  of  flame  thickness,  meaning  the  bottom  curve  in  Figure  3.25  represents  the  minimum 
flame  thickness  achievable  for  this  parameter  set.  Figure  3.26  shows  how  Xi  varies  along  this  curve  as  a 
function  of  Ea.  The  flame  thickness  decreases  with  increasing  activation  energy,  but  not  quickly  enough  to 
fall  close  to  the  experimental  value  of  approximately  Xi  =  0.25  cm.  The  laminar  flames  computed  using  this 
approach  are  approximately  2-3  times  too  thick  within  the  range  of  reasonable  activation  energies  for  this 
mixture.  The  choice  of  Ea,  is  then  somewhat  arbitrary,  so  we  set  Ea  =  67.551?To,  the  same  value  used  for 
the  stoichiometric  mixture,  and  choose  n0  =  4.71  x  10-5  g/s-cm-K0  7  and  A  =  1.053  x  1012  cm3/g/s,  which 
gives  the  computed  Si  =  15.7  cm/s  and  Xd  =  0.69  cm.  These  parameters  and  the  corresponding  computed 
flame  and  detonation  properties  are  given  in  the  second  column  of  Table  6. 

We  then  use  both  models  for  13%  methane-air  combustion  to  compute  flame  acceleration  in  configuration 
17.4  with  £  =  0.3.  The  resulting  flame  velocities  are  compared  with  velocities  measured  in  experiments  with 
12%  and  13%  methane  concentrations  [72]  in  Figure  3.27.  We  find  a  significant  difference  in  the  flame 
acceleration  between  the  two  cases.  This  observation  is  in  spite  of  the  fact  that  both  models  produce  similar 
Si  and  T/.  The  principle  difference  between  the  two  cases  is  the  thickness  of  the  computed  flames.  The  model 
2  flame  is  10  times  thicker  than  the  model  1  flame.  The  result  is  a  much  more  slowly  accelerating  turbulent 
flame  brush.  The  velocities  computed  for  the  model  2  flame  do,  however,  agree  reasonably  well  with  the 
experimental  measurements.  Neither  model  predicts  DDT,  nor  was  it  observed  in  experiments  either  for  this 
rich  case. 

3.3.3.  Discussion  and  Other  Approaches  Under  Consideration 

Figure  3.27  shows  that  there  is  some  sensitivity  in  flame  acceleration  to  the  thickness  of  the  modeled  flame. 
This  is  in  contrast  to  the  lack  of  any  pronounced  sensitivity  to  other  flame  parameters,  such  as  Si  and  Tj,. 
Further  testing  is  needed  to  isolate  this  effect  by  considering  a  series  of  flames  with  the  same  Si  and  T/,  but 
different  Xi,  controlled  by  varying  k.  Physically,  there  may  be  some  justification  for  rich  mixtures  having 
larger  k.  The  thermal  diffusivity  of  methane  is  slightly  larger  than  that  of  air,  so  the  additional  methane 
present  in  rich  mixtures  and  product  gases  leads  to  larger  effective  diffusivities.  The  calibrated  k  for  model 
2  is  still  much  larger  than  one  would  expect,  however. 
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Table  6.  Material,  chemical,  and  reaction-wave  parameters  for  rich  (13%)  methane-air 


Quantity 

Value 

Value 

Definition 

Input 

Po 

1  atm 

1  atm 

Initial  pressure 

T0 

298  K 

298  IC 

Initial  temperature 

Po 

1.1042  x  10~3  g/cm3 

1.1042  x  10~3  g/cm3 

Initial  density 

7 

1.19 

1.27 

Adiabatic  index 

Ms 

27 

27 

Molecular  weight 

A 

2.6  x  1012  cm3/g/s 

1.053  x  1012  cm3/g/s 

Pre-exponential  factor 

E„ 

57.55  7?  To 

67.55  7?  T0 

Activation  energy 

q 

36.3  R  T0  /M 

26.71 7?T0/M 

Chemical  energy  release 

vo 

3.6  x  10~6  (g/s-cm-K0-7) 

3.6  x  10-6  (g/s-cm-K0  7) 

Transport  constants 

«o  =  Do 

6.25  x  10-6  (g/s-cm-K0-7) 

4.71  x  10"5  (g/s-cm-K0  7) 

Transport  constants 

Output 

Si 

24.1  cm/s 

15.7  cm/s 

Laminar  flame  speed 

Tb 

2025  K 

1990  K 

Adiabatic  flame  temperature 

Xl 

0.0631  cm 

0.698  cm 

Laminar  flame  thickness 

Dcj 

1728  m/s 

1800  m/s 

C-J  detonation  velocity 

Xd 

0.54  cm 

0.69  cm 

ID  half-reaction 

detonation  thickness 

The  single  component  reaction  model,  in  which  fuel  is  converted  completely  into  product  in  a  single 
reaction  step,  is  not  correct  for  rich  mixtures  in  which  there  is  fuel  left  over  after  oxidation  in  complete, 
i.e. ,  when  the  oxidizer  is  depleted.  For  this  situation,  the  oxidizer  is  the  limiting  component  and  should  be 
explicitly  computed.  Such  a  two-component  reaction  model  and  its  effect  on  flame  acceleration  and  DDT 
will  be  discussed  in  more  detail  in  the  following  section. 
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Figure  3.21.  (a)  Dcj  and  Tf,  computed  using  chemical  equilibrium  software  [73,  33]  and  (b)  A“  and  Si 
measured  experimentally  [14,72]  as  a  function  of  methane  concentration  for  laminar  flames  and  detonations. 


(a)  (b) 
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Figure  3.22.  Parametric  curves  for  which  (a)  Tf,  =  1971  K,  Dqj  =  1692  nr/s  and  (b)  S/  =  23.7  cm/s, 
Xd  =  0.8,  and  1.6  cm  for  a  lean  (7%)  methane-air  mixture  computed  using  the  single-component  reaction 
model  (Eq.  (33)).  The  point  of  intersection  in  (a)  gives  the  values  of  qM/RT0  and  7  used  to  compute  the 
curves  in  (b).  Two  curves  for  Xd  are  given  to  indicate  uncertainty  in  correlation  between  detonation  cell  size 
and  Xd- 
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Figure  3.23.  Velocity  of  the  leading  edge  of  the  calculated  7%  methane-air  flame  (<E>(/>  =  0.72)  as  a  function 
of  its  location  in  the  channel  in  Configuration  17.4,  £  =  0.3  for  (a)  various  values  of  Si  and  fixed  T),  and 
(b)  various  values  of  Tj  and  fixed  Si  along  with  experimental  measurements  for  6.5%  and  8%  methane-air 
flames  in  a  similar  configuration  [72]. 


Figure  3.24.  Parametric  curves  for  which  (a)  X),  =  1990  K,  Dcj  =  1800  m/s  and  (b)  Si  =  23.7  and 
15.7  cm/s,  Xd  =  0.69  and  1.1  cm  for  a  rich  (13%)  methane-air  mixture  computed  using  the  single-component 
reaction  model  (Eq.  (33)).  The  point  of  intersection  in  (a)  gives  the  values  of  qM/ RTq  and  7  used  to  compute 
the  curves  in  (b).  Two  curves  for  Xd  are  given  to  indicate  uncertainty  in  correlation  between  detonation  cell 
size  and  Xd- 
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Figure  3.25.  Variation  of  k0  needed  to  attain  Si  =  15.7  or  23.7  cm/s  using  the  values  of  A  for  which 
Xd  =  0.69  or  1.1  cm  shown  in  Figure  3.23b  as  a  function  of  Ea/RTo- 


Figure  3.26.  Laminar  flame  thickness  computed  using  values  of  k0  on  the  lowest  curve  in  Figure  3.24 
(Si  =  15.7  cm/s,  Xd  =  0.69  cm)  as  a  function  of  Ea/RT0. 


67 


Figure  3.27.  Velocity  of  the  leading  edge  of  the  calculated  13%  (</>  =  1.43)  methane-air  flame  as  a  function 
of  its  location  in  the  channel  in  Configuration  17.4,  £  =  0.3  computed  using  the  two  parameter  sets  given 
in  Table  6  along  with  experimental  measurements  for  12%  and  13%  methane-air  mixtures  in  a  similar 
configuration  [72]. 
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3.4.  Propagation  in  a  Gradient  of  Fuel  Composition  (Task  3) 

While  DDT  and  detonation  propagation  in  uniform  fuel-air  mixtures  has  been  studied  extensively,  much  less 
is  known  about  the  behavior  of  deflagration  waves,  DDT,  and  detonations  in  systems  where  the  fuel  and  air 
are  only  partially  premixed,  that  is,  where  the  local  fuel  concentration  is  not  uniform.  Such  heterogeneous 
methane-air  systems  may  form  in  the  unventilated,  sealed  areas  of  underground  coal  mines.  If  methane  gas 
seeps  in  from  the  tunnel  roof,  then  a  gradient  will  exist  from  roof  to  floor  while  the  methane  diffuses  into  the 
air.  The  magnitude  of  this  gradient  will  depend  on  the  methane  influx  rate  from  the  tunnel  roof.  If  air  leaks 
through  seals  into  a  methane-rich  sealed  area,  then  buoyancy  will  cause  the  inflowing  air  to  collect  near  the 
tunnel  floor  and  the  existing  methane-rich  atmosphere  to  collect  near  the  tunnel  roof  until  diffusion  leads 
to  their  complete  mixing.  In  this  case,  the  magnitude  of  the  composition  gradient  will  depend  on  the  air 
inflow  rate  into  the  sealed  area.  Understanding  how  a  flame  accelerates,  whether  it  can  undergo  transition 
to  a  detonation,  and  how  that  detonation  propagates  are  important  general  problems  anywhere  natural  gas 
accumulates. 

In  premixed  systems,  the  propagation  speed  and  the  cellular  structures  formed  by  instability  of  a  det¬ 
onation  front  depend  strongly  on  the  equivalence  ratio  of  the  fuel-air  mixture.  The  general  trend  is  for  the 
detonation  to  slow  and  cell  sizes  to  grow  as  the  mixture  moves  away  from  stoichiometric.  In  methane-air 
systems,  the  detonation  cell  width  was  found  to  vary  between  21  and  51  cm  over  a  methane  concentration 
range  of  only  3.5%  [72].  These  observations  suggest  that  in  partially-premixed  systems  where  the  concentra¬ 
tion  of  fuel  is  nonuniform,  both  the  propagation  speed  and  the  detonation  cell  size  can  vary  as  a  function  of 
space.  The  implications  of  this  spatial  dependence  have  been  studied  in  systems  where  the  gradient  of  the 
fuel  concentration  is  parallel  to  the  direction  of  propagation  [86-88].  The  magnitude  of  the  concentration 
gradient  was  found  to  be  a  determining  factor  in  whether  or  not  a  detonation  could  propagate  through  the 
nonuniform  mixture. 

Less  work  has  been  done  on  systems  for  which  the  gradient  in  fuel  concentration  is  normal  to  the 
direction  of  propagation  of  the  detonation.  Experiments  have  shown  that  the  sizes  of  detonation  cells  vary 
spatially  in  relation  to  the  variation  in  fuel  concentration  [89]  and  that  the  shape  of  the  detonation  wave 
is  curved  in  these  systems  [89,  90].  Calhoon  and  Sinha  [91]  calculated  the  structure  of  a  detonation  in  the 
mixing  layer  of  laminar,  initially  non-premixed,  co-flowing  streams  of  fuel  and  air.  They  found  a  complex 
wave  structure  that  resembled  a  typical  triple  flame  found  in  low-speed  mixing  layers  [92,  93].  A  detonation 
formed  near  the  line  of  stoichiometry  in  the  mixing  layer  that  extended  into  the  fuel-rich  and  fuel-lean 
mixtures  above  and  below  this  point.  In  these  regions,  the  detonation  was  curved  due  to  variations  in  the 
local  propagation  speed.  Far  from  the  line  of  stoichiometry  where  the  mixture  was  either  too  rich  or  too 
lean  to  support  a  detonation,  the  leading  shock  detached  from  the  reaction  zone  leaving  behind  high-speed 
deflagrations.  Downstream,  a  diffusion  flame  consumed  the  excess  fuel  and  oxidizer  that  were  not  burned  in 
the  fuel-rich  and  fuel-lean  detonations  and  deflagrations. 

Here  we  generalize  this  to  a  partially  premixed  system  (similar  to  that  used  by  Calhoon  and  Sinha  [91]) 
in  which  the  local  methane-air  equivalence  ratio  varies  as  a  function  of  space  and  time.  Ultimately,  we  wish 
to  consider  flame  acceleration  and  DDT  in  configurations  similar  to  those  considered  in  Sections  3. 1-3.3. 
This  work,  however,  is  still  ongoing,  and  in  the  following  sections  we  give  only  an  overview  of  the  calibration 
of  the  numerical  and  physical  models  being  used  .  We  then  test  these  models  by  performing  calculations  of 
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a  detonation  propagating  in  a  heterogeneous  methane-air  mixture.  These  preliminary  results  are  discussed 
in  context  of  the  previous  experimental  and  numerical  studies  described  above.  Eventually,  we  plan  to  use 
the  same  models  to  compute  flame  acceleration  and  transition  to  detonations  in  obstructed  channels. 


3.4.1  The  Two-Component  Model 

Since  now  the  relative  concentrations  of  fuel  and  oxidizer  can  vary  throughout  the  system,  we  must  consider 
a  two-component ,  single-step  reaction  model,  vpF  +  voO  — >  (vp)P,  in  which  vp  moles  of  fuel  ( F )  combine 
with  vo  moles  of  oxidizer  ( O )  to  form  vp  moles  of  product  (P).  The  molecular  weights  of  the  fuel  and 
oxidizer  species  are  given  by  Wf  and  Wo,  respectively.  In  order  to  compute  partially  premixed  flames  and 
detonations,  it  is  necessary  to  consider  the  dynamics  of  both  the  fuel  and  oxidizer  species,  the  mass  fractions 
of  which  (Y  and  X ,  respectively)  are  governed  by 


d(pY) 

-W^  +  X7  ■(PYv)~X7  ■(pDfX7Y)-Pn  =  0,  (3.11) 

d(pX)  +  X7  ■  (pXv)  -  X  ■  (pD0\7X)  ~  pSCl  =  0  .  (3.12) 

at 

These  equations  are  solved  concurrently  with  Eqs.  (2. 1-2.5).  Here,  we  assume  a  two-component  reaction 
model  for  which  we  write  the  fuel  consumption  rate, 


Cl  =  -pA{(l))YXexp(-Ea(<j))/RT),  (3.13) 

that  explicitly  depends  on  the  local  concentrations  of  both  fuel  and  oxidizer.  The  pre-exponential  factor  A 
and  the  activation  energy  Ea  are  assumed  to  vary  with  the  local  equivalence  ratio  4>-  In  order  to  ensure  the 
model  gives  the  correct  flame  temperature  and  flame  speed  as  a  function  of  equivalence  ratio,  we  also  allow 
q  and  k0  to  vary  as  functions  of  <fi .  For  a  two-component  mixture,  we  define  a  local  mixture  fraction  Z  that 
does  not  vary  across  the  reaction  zone  [94,  95], 

=  SY-X-  SYlean  +  1 

1  +  S-Xrich-SYlean'  1  ‘  ’ 

where  S  =  sY^/X^  is  the  global  equivalence  ratio,  s  =  v0W0/vfWf  is  the  stoichiometric  coefficient  of  the 
reaction,  Yiean  is  the  mass  fraction  of  fuel  in  the  leanest  region  of  the  flow  (where  X  =  Xca)1  and  Xrich  is 
the  oxidizer  mass  fraction  in  the  richest  region  of  the  flow  (where  Y  =  Yoo)-  Since  Z  is  a  conserved  quantity 
across  the  flame  front,  it  can  be  related  to  the  fuel  and  oxidizer  mass  fractions  in  the  unburned  mixture  just 
upstream  of  the  reaction  zone  according  to 


Yunb  —  Z(1  —  Yiean )  +  Y( 


ear n 


Xuni,  —  (1  Z)(l  Xrich)  T  Xrich. 


(3.15) 


The  local  equivalence  ratio  can  then  be  computed  based  on  these  unburned  mass  fractions, 

0 _  S'  Yu  f,}, 

Xunb 


(3.16) 


We  define  X ^  and  Y^  as  the  maximum  possible  oxidizer  and  fuel  mass  fractions  in  the  domain,  re¬ 
spectively.  In  the  equations  above,  the  fuel  and  oxidizer  mass  fractions,  Y  and  X,  have  been  normalized  by 
these  values  so  that  they  vary  between  0  and  1. 
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3.4.1a.  Calibration  for  Stoichiometric  Mixture 


Table  7.  Input  model  parameters  and  computed  properties  of  reaction  waves  for  five  different 
methane-air  mixtures  using  a  two-component  reaction  model. 


Quantity 

5.9% 

(</>  =  0.6) 

7% 

((t>  =  0.72) 

9.5% 

(0=1) 

11.2% 

(0=1.2) 

13% 

(0  =  1-43) 

Input 

P0  (atm) 

1 

1 

1 

1 

1 

r0(K) 

298 

298 

298 

298 

298 

Po  (g/cm3) 

1.1042  x  10"3 

1.1042  x  10"3 

1.1042  x  10~3 

1.1042  x  10~3 

1.1042  x  10"3 

7 

1.172 

1.18 

1.197 

1.244 

1.27 

Ms  (g/rnol) 

27 

27 

27 

27 

27 

A  (cm3/g-s) 

1.0697  x  1011 

7.0x10“ 

1.78xl013 

7.583  x  1014 

9.343xl012 

Eat  RTo 

42.0 

52.55 

67.55 

95 

77.55 

qM/RTo 

33.66 

36.81 

39.0 

37.2 

38.29 

vq  (g/s-cm-K0-7) 

3.6xl0"6 

3.6xl0"6 

3.6xl0-6 

3.6xl0"6 

3.6xl0~6 

Ko  =  Df0  =  Do0 
(g/s-cm-K07) 

1.53  x  10"5 

2.1xl0"5 

3.25xl0-5 

5.92  x  10"5 

5.25xl0"5 

Output 

Si  (cm/s) 

15.4 

21.6 

40.3 

38.0 

15.9 

Tb  (I<) 

1770 

1971 

2210 

2110 

1990 

xi  (cm) 

0.24 

0.24 

0.21 

0.37 

0.76 

DCj  (m/s) 

1587 

1691 

1820 

1827 

1800 

xd  (cm) 

2.02 

1.32 

0.23 

0.29 

0.39 

Calibration  of  the  model  parameters  is  performed  in  a  manner  similar  to  that  described  in  Sections  3.1.1  and 
3.3.2.  Target  laminar  flame  properties  (flame  speed,  adiabatic  flame  temperature,  laminar  flame  thickness) 
and  detonation  properties  (Chapman-Jouget  velocity,  ID  detonation  thickness)  are  given  in  Table  1  for 
several  different  metlrane-air  mixtures.  We  solve  the  steady-state  laminar  flame  equations  (Eqs.  (3.5),  (3.6)) 
and  the  one-dimensional  ZND  equations  for  a  propagating  detonation  wave  (Eqs.  (3.7)-(3.10))  using  the 
two-component  reaction  model  (Eq.  3.13)  to  obtain  a  set  of  parameters  for  which  Tb,  Si,  Xi,  Dqj,  and  Xd 
fall  within  these  target  ranges  for  each  of  the  three  compositions.  Since  Tf,  and  Dqj  do  not  depend  on  the 
reaction  model,  the  values  for  7  and  q  are  chosen  as  the  same  points  of  intersections  of  the  curves  shown  in 
Figures  3.1a,  3.22a,  and  3.24a.  These  values  are  listed  in  Table  7.  For  this  two-component  reaction  model, 
the  actual  magnitude  of  the  non-dimensional  heat  release  q*  is  scaled  by  so  that  q  =  q*Y00  and  the 
magnitude  of  the  actual  pre-exponential  factor  A*  is  scaled  by  X x  so  that  A  =  A*X00. 

Next,  we  next  calibrate  A  and  Ea  for  each  composition  based  on  the  Si  and  Xd  given  in  the  table. 
The  use  of  the  two-component  reaction  model  causes  only  small  differences  in  the  computed  detonation 
structure  from  that  computed  using  the  single-component  model.  The  choice  of  reaction  model,  however, 
significantly  affects  the  speed  and  thickness  of  a  laminar  flame.  Values  of  A  and  Ea  that  give  Si  and  Xd 
in  the  target  range  for  a  stoichiometric  mixture  are  plotted  in  Figure  3.28.  These  curves  do  not  intersect 
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within  the  reasonable  range  of  Ea.  We  can  alter  the  flame  speed  for  fixed  A ,  Ea  by  changing  the  thermal 
conductivity  of  the  mixture.  The  approach  we  take  here  is  to  choose  a  point  in  the  A-Ea  parameter  space 
that  gives  Xd  within  the  target  range  shown  in  Table  1.  In  principle,  we  could  choose  any  such  A,  Ea  pair, 
but  to  make  direct  comparisons  with  results  obtained  using  the  single-component  reaction  model  we  choose 
the  same  Ea/RT0  =  67.55  and  select  A  =  1.78  x  1013  cm3/g-s  such  that  the  computed  Xd  =  0.229  cm.  This 
detonation  thickness  is  identical  to  that  obtained  using  the  single-component  reaction  model  as  shown  in 
Figure  3.29.  Given  this  set  of  parameters,  g,  7,  A,  Ea,  we  then  increase  k0  to  give  Si  in  the  target  range. 
The  resulting  na  =  3.25  x  10-5  g/s-cm-K0-'  is  much  larger  than  the  actual  conductivity  of  the  methane-air 
mixture,  but  when  used  along  with  the  rest  of  the  parameter  set,  it  gives  the  correct  laminar  flame  speed 
and  a  laminar  flame  thickness  that  is  comparable  to  experimental  measurements.  Comparisons  between  two- 
component  model  flame  structures  and  the  single-component  model  flame  structure  are  made  in  Figure  3.30. 
The  flames  represented  by  the  blue  and  black  lines  both  propagate  at  the  same  speed,  but  are  quite  different 
in  terms  of  thickness. 

We  next  compare  flame  acceleration  computed  using  single-component  and  two-component  reaction 
models.  Two-dimensional  simulations  of  configuration  17.4  (cf.  Section  3.1)  are  performed  for  a  blockage 
ratio  £  =  0.3.  Grid  sizes  are  chosen  so  that  both  the  smallest  dxmin  and  largest  dxmax  are  equal  to  those 
shown  in  Section  3.2  to  produce  converged  results,  i.e.,  dxmin  =  0.009062  cm  and  dxmax  =  0.145  cm.  Flame 
velocities  as  a  function  of  position  in  the  channel  are  shown  for  three  different  cases  in  Figure  3.31.  The 
single-component  case  with  parameters  shown  in  Table  2  is  given  by  the  red  line.  The  laminar  flame  speeds 
are  identical  for  the  two  cases  shown  for  the  two-component  model  (green  and  blue  lines)  ,  but  the  case 
with  larger  k0  and  smaller  A  (blue  line)  gives  a  much  larger  flame  thickness  (see  Figure  3.30.)  This  leads  to 
a  restriction  in  the  rate  of  growth  of  total  flame  surface  area  due  to  wrinkling,  and  hence  a  slower  energy 
output  and  flame  acceleration.  This  curve  seems  to  more  closely  match  the  experimental  data  shown  in  the 
figure  (light  blue,  purple)  than  either  curve  generated  with  smaller  k0- 

3.4.1b.  Calibration  for  Non-Stoichiometric  Mixtures 

The  calibration  of  non-stoichiometric  values  of  </>  can  be  done  in  a  variety  of  ways.  Here,  we  discuss  three 
different  approaches.  First,  we  assume  that  the  transport  and  reaction  parameters  are  constant  over  the 
entire  range  of  (f>.  Second,  we  perform  the  same  calibration  procedure  used  above  to  compute  the  parameters 
for  a  stoichiometric  mixture  to  find  sets  of  parameters  for  two  additional  equivalence  ratios  and  then  use  a 
second-order  polynomial  to  interpolate  the  results  to  other  values  of  4>.  Third,  we  calibrate  sets  of  parameters 
for  four  additional  equivalence  ratios  (for  a  total  number  of  calibration  mixtures  of  five)  and  use  fourth-order 
polynomials  to  construct  the  variation  of  the  reaction  parameters  with  cf>  for  all  intermediate  values.  Each 
of  these  approaches  will  be  discussed  in  more  detail  below. 

The  most  straightforward  approach  is  to  use  constant  values  of  the  reaction  parameters  g*,  A*,  7,  Ea, 
and  k0.  The  terms  q  =  q*Yoc  and  A  =  A*^^  that  appear  in  Eqs.  (3.11),  (3.12),  and  (2.3)  naturally  vary 
with  (f)  and  thus  allow  for  the  detonation  parameters  to  change  with  the  mixture  composition.  The  values 
of  these  parameters  calibrated  for  a  stoichiometric  mixture  are  used  for  a  range  of  cf>  from  0.6-1. 4,  and 
Figures  3.33a  and  b  show  the  resulting  variations  of  Dqj  and  Xd  computed  using  this  approach.  These 
results  are  compared  with  accepted  values  of  Dqj  obtained  using  a  chemical  equilibrium  code  [73]  and  of  Xd 
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based  on  the  correlation  [59]  between  reaction  zone  width  and  measured  detonation  cell  sizes  [14,  72]  used 
above.  The  computed  values  of  Dqj  and  Xj  deviate  substantially  from  the  accepted  values  for  0  ^  1. 

In  the  second  approach,  we  attempt  to  improve  the  predictive  capability  of  the  reaction  model  for  lean 
and  rich  mixtures  by  allowing  the  reaction  parameters  to  vary  with  0.  Two  additional  values  of  0  are 
chosen  for  calibration,  one  lean  (0  =  0.72,  7%  methane)  and  one  rich  (0  =  1.43,  13%  methane),  with  target 
detonation  and  flame  properties  given  in  Table  1.  The  procedure  used  to  find  A,  Ea,  q,  7,  and  kq  for  these 
two  equivalence  ratios  is  the  same  as  that  used  to  calibrate  the  stoichiometric  mixture.  These  parameters 
are  also  summarized  in  Table  7.  In  performing  these  calibrations,  we  adjusted  the  parameters  so  that  the 
computed  values  of  Xi  varied  with  equivalence  ratio  in  a  manner  similar  to  that  observed  in  experiments. 
The  value  computed  for  0  =  1  is  somewhat  larger  than  the  target  value  shown  in  Table  1.  Values  of  Ea, 
«o,  and  A  for  the  7%  and  13%  calibration  mixtures  were  adjusted  to  ensure  both  the  lean  and  rich  flames 
were  thicker  than  the  stoichiometric  flame,  in  accordance  with  experimental  measurements.  The  parameter 
sets  computed  and  given  in  Table  7  do  not  constitute  a  definitive  reaction  model  for  the  respective  mixtures 
since  there  is  some  uncertainty  in  the  experimental  measurements  used  to  calibrate  them.  Hence,  a  variety 
of  possible  model  parameter  sets  exist  that  could  adequately  reproduce  the  data  given  in  Table  1. 

We  then  fit  a  second  order  polynomial  to  the  calibrated  values  of  A,  Ea,  g,  Kq,  and  7  to  give  a  continuous 


set  of  values  for  these  parameters  as  a  function  of  0, 

A  =  Ast  exp  [— 18.38802  +  43.1850  -  24.796] ,  (3.17) 

Ea  =  K,st(-O.632102  +  1.88030  -  0.2482),  (3.18) 

q  =  9st(— O.342702  +  0.79050  +  0.5522),  (3.19) 

7  =  75T(O.128302  -  0.170  +  1.0417),  (3.20) 

Ho  =  kg,st(  1.37070  -  0.352),  (3.21) 


where  the  subscript  st  indicates  the  parameter’s  value  in  the  model  stoichiometric  mixture  given  in  Table  7. 

This  particular  choice  of  fitting  curve  is  empirical  and  was  made  to  most  closely  match  the  data  for 
Dqj  and  xj-  The  transformation  of  A  to  logarithmic  coordinates  was  done  to  ensure  that  A  cannot  become 
negative  for  any  value  of  0.  Values  of  A,  Ea,  g,  Ko,  and  7  calculated  using  Eqs.  (3.17-3.21)  are  shown  in 
Figures  3.32a-e  for  a  range  of  0,  and  the  values  of  Dqj  and  xj  computed  using  these  parameters  are  shown 
in  Figures  3.33a  and  b.  This  second  order  approach  does  a  better  job  at  describing  both  Dqj  and  Xj  for 
lean  values  of  0,  but  it  introduces  a  non-physical  local  minimum  of  Dqj  for  rich  mixtures.  Additionally,  the 
maximum  Dqj  occurs  for  0  =  1  in  contrast  to  the  exact  Dqj,  which  peaks  between  0  =  1.1  and  1.2. 

The  third  approach  addresses  this  shortcoming  by  adding  more  control  points,  i.e,  calibrated  parameter 
sets  for  more  mixture  compositions,  and  seeking  a  higher-order  polynomial  fit  for  the  reaction  parameters. 
Two  additional  values  of  0  are  considered:  0  =  0.6  and  0  =  1.2.  The  same  calibration  process  used  to  obtain 
the  transport  and  reaction  coefficients  for  0  =  0.72,  1,  and  1.43  is  used  to  calibrate  0  =  0.6,  and  0  =  1.2. 
We  found  it  impossible  to  obtain  a  set  of  parameters  with  realistic  values  of  activation  energy  that  produces 
Xd,  Si,  and  Xi  that  are  all  close  to  the  accepted  values.  Instead,  we  adjusted  A,  Ea,  and  nQ  to  match  xj  and 
Si  and  maintain  the  correct  qualitative  behavior  of  Xi  in  relation  to  the  values  obtained  for  the  other  control 
points.  The  transport  and  reaction  coefficients  obtained  for  0  =  0.6  and  1.2  are  shown  in  Table  7. 
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We  then  fit  fourth-order  polynomials  to  all  five  of  the  control  points  for  the  various  reaction  parameters, 


resulting  in  the  following  set  of  equations: 

A  =  Ast  exp  [-225.62O104  +  836.25640s  -  1131.262102  +  674.69150  -  154.0657] ,  (3.22) 

Ea  =  £0,st(-24.56504  +  92. 8880s  -  127.5902  +  76.4840  -  16.216),  (3.23) 

q  =  gSr(2.358504  -  7.8387 03  +  8.445902  -  2.93940  +  0.9737),  (3.24) 

7  =  7st(-1.55O904  +  5.985403  -  8.2683 02  +  4.92470  -  0.091),  (3.25) 

k0  =  kg,st(- 42.7O704  +  161.O403  -  219.O702  +  129.130  -  27.391),  (3.26) 


which  are  plotted  along  with  the  second-order  fits  in  Figures  3.32a-e.  The  additional  control  points  add  local 
maxima  and  minima  to  the  curves  for  the  individual  reaction  parameters,  which  have  no  reasonable  physical 
description.  They  appear  as  a  result  of  trying  to  force  a  single  reaction  model  to  describe  two  disparate 
burning  regimes:  low-speed  flames  and  high-speed  detonations.  Taken  collectively  and  used  in  conjunction 
with  Eq.  (3.13),  these  parameters  produce  Dqj  and  Xd  shown  in  Figures  3.33a  and  b  that  exhibit  the  desired 
features  of  the  exact  Dcj  and  Xd ■  Similarly  good  agreement  is  obtained  for  Tj,  and  Si  (not  shown),  and 
the  variation  in  xi  is  qualitatively  correct  although  each  individual  value  is  larger  than  the  exact  value.  We 
caution,  however,  that  Eqs.  (3.22-3.26)  should  not  be  used  outside  of  the  range  of  0  =  0.6-1.43  without 
further  calibration. 

3.4.2.  Detonations  in  Uniform,  High-Activation-Energy  Mixtures 

The  ultimate  goal  of  this  work  is  to  use  the  model  developed  in  Section  3.4.1  to  compute  turbulent  flame 
acceleration,  shock-flame  interaction,  DDT,  and  detonation  propagation  in  partially-premixed  systems  in 
which  a  bulk  gradient  in  the  local  fuel  concentration  exists.  This  work  is  ongoing,  and  here  we  only  present 
preliminary  results  concerning  the  propagation  of  a  detonation  through  such  a  concentration  gradient. 

Consider  a  nonuniform  mixture  of  methane  and  air  in  a  32  cm-wide  by  1024  cm-long  channel  as  shown 
in  Figure  3.34.  The  mixture  at  the  top  of  the  channel  (y  =  32  cm)  is  taken  to  be  fuel  rich,  0r,  and  that 
at  the  bottom  of  the  channel  (y  =  0  cm)  is  taken  to  be  fuel  lean,  0;.  The  initial  concentration  gradient  is 
set  by  assuming  a  linear  variation  in  0  from  0;  to  0r.  The  mixture  is  ignited  directly  by  a  strong  shock 
with  Mach  number,  Ms,  equal  to  10.  The  shock  is  placed  initially  2  cm  to  the  right  of  the  closed  end  of 
the  channel.  The  temperature  of  the  compressed  gas  in  the  region  behind  the  shock  is  sufficiently  high  to 
allow  a  detonation  to  ignite  in  a  short  period  of  time,  which  then  propagates  from  left  to  right  through  the 
domain. 

We  consider  first  the  propagation  of  detonation  waves  through  several  different  uniform  methane-air 
mixtures,  ranging  in  equivalence  ratio  from  0  =  0.7  to  0  =  1.1.  For  all  cases,  we  set  0j  =  0r  so  that  there  is 
no  concentration  gradient  within  the  mixtures.  The  detonation  propagation  velocities  in  each  mixture  are 
shown  in  Figure  3.35  as  a  function  of  the  x-location  of  the  leading  edge  of  the  reaction  zone  in  the  channel. 
The  thinner  horizontal  lines  indicate  the  value  of  Dqj  for  each  mixture.  After  an  initial  adjustment  period, 
equal  to  about  15  channel  widths,  the  propagation  velocity  decays  to  a  level  consistent  with  Dcj-  There  is 
little  difference  between  the  propagation  velocities  for  the  0  =  1  and  1.1  cases  shown  in  the  figure,  consistent 
with  the  behavior  of  Dcj  over  that  range  of  0,  but  there  is  a  noticeable  drop  in  velocity  for  the  lean  cases. 
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The  leanest  case,  </>  =  0.7,  cannot  support  a  detonation  in  this  channel,  and  a  steep  drop  in  flame  velocity 
occurs  when  the  the  leading  shock  and  flame  decouple  after  propagating  together  for  about  180  cm  (or  about 
6  channel  widths).  This  is  not  surprising  since  the  detonation  cell  size  for  a  mixture  with  <f>  =  0.7  is  close 
to  90  cm  [14,  72],  which  is  much  larger  than  the  size  of  the  channel.  The  observed  detonation  failure  is  not 
related  to  any  thermal  or  mechanical  losses  that  are  not  included  in  our  model.  The  failure  occurs  because 
of  the  strong  instabilities  of  the  detonation  front  in  a  system  with  high  activation  energy.  Both  transverse 
and  longitudinal  instabilities  contribute  to  the  spatial  and  temporal  fluctuations  of  the  detonation  velocity. 
When  the  front  slows  down  to  Dcj,  and  only  a  few  triple  points  remain  in  the  channel,  these  points  can 
easily  weaken  as  a  result  of  a  fluctuation.  In  a  system  with  high  activation  energy,  the  induction  delay  behind 
a  shock  rapidly  decreases  as  the  shock  weakens.  Thus,  the  reaction  in  the  vicinity  of  weakened  triple  points 
may  quench,  and  the  flame  can  decouple  from  the  leading  shock  as  a  result  of  a  fluctuation.  In  a  very  long 
channel  without  losses,  the  detonation  would  eventually  reignite,  but  this  may  never  happen  in  channels 
limited  to  practical  lengths. 

Figure  3.36a  shows  a  sequence  of  density  plots  in  the  vicinity  of  the  propagating  reaction  front  as  the 
detonation  transitions  to  a  shock-flame  complex.  The  red  parts  are  high-density  regions  containing  unreacted 
material,  which  was  compressed  by  an  inert  shock  wave  that  decoupled  from  the  reaction  front.  Green  parts 
of  the  front  are  overdriven  detonations.  Several  transverse  instabilities  on  the  detonation  front  are  visible 
prior  to  the  detonation  failure.  These  triple  points  are  the  regions  of  high  pressure  that  produce  the  cellular 
patterns  observed  in  experiments.  Numerically,  we  track  the  trajectories  of  these  triple  points  by  recording 
the  maximum  pressure  experienced  in  each  computational  cell.  An  example  of  such  a  numerical  smoke  foil 
for  the  uniform  <f>  =  0.7  detonation  is  shown  in  Figure  3.37a.  The  highest  pressures  are  represented  by 
the  black  traces  in  the  figure.  Near  the  left  side  of  the  smoke  foil,  as  the  detonation  is  just  beginning  to 
propagate,  small-scale  instabilities  develop  along  the  initially  planar  and  smooth  detonation  front.  These 
instabilities,  most  likely  triggered  by  numerical  noise  in  the  integration  routine,  form  cells  that  grow  from 
small  and  regular  diamonds  when  the  detonation  is  initially  overdriven  to  large,  irregular  quadrilaterals  as 
the  detonation  approaches  the  CJ  velocity.  In  the  vicinity  of  180  cm  from  the  closed  end  of  the  channel  the 
triple  points  weaken  and  disappear  as  the  detonation  fails. 

A  similar  process  occurs  for  the  <j>  =  0.8  case,  but  the  detonation  is  able  to  survive  and  propagate 
at  near-CJ  speed  for  a  much  longer  distance  than  the  detonation  in  the  <j>  =  0.7  mixture.  Shock-flame 
decoupling  occurs  around  320  cm  from  the  closed  end  wall.  Again,  triple  points  all  weaken  and  die  out  as 
they  propagate  along  the  detonation  front  (Figure  3.37b).  This  is  also  to  be  expected  since  the  cell  size  for 
a  detonation  in  a  </>  =  0.8  mixture  is  approximately  65  cm  [14,  72],  which  is  also  larger  than  the  height  of 
the  channel. 

For  the  <j>  =  0.9  and  1.0  cases,  the  detonations  are  able  to  survive  and  propagate  indefinitely  through 
the  channel  at  average  speeds  roughly  equal  to  the  CJ  velocities  for  each  mixture.  It  is  difficult  to  get 
a  precise  measurement  of  the  propagation  velocity  because  of  the  influence  of  longitudinal  and  transverse 
instabilities.  For  high-activation  energy  mixtures,  the  longitudinal  propagation  velocity  has  been  shown 
to  fluctuate  about  the  CJ  velocity,  locally  weakening  and  strengthening  the  detonation.  The  transverse 
instabilities  that  give  rise  to  triple  points  also  cause  ambiguity.  At  some  times  and  locations,  the  leading 
edge  of  the  front  between  two  triple  points  is  an  overdriven  detonation  with  velocity  greater  than  Dqj,  and 
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at  other  times  and  locations  a  decoupled  shock- flame  complex  with  velocity  less  than  Dqj-  A  sequence 
of  density  fields  taken  from  the  </>  =  1.0  calculations  at  several  instants  over  a  characteristic  detonation¬ 
cell  crossing  time  are  shown  in  Figure  3.36b.  The  color  mappings  have  the  same  interpretation  as  those 
given  above  for  the  <f>  =  0.7  case.  While  there  are  intermittent  periods  of  shock- flame  decoupling  along  the 
reaction  front,  there  is  never  a  global  shock-flame  decoupling  like  that  shown  in  Figure  3.36a.  Instead,  a 
myriad  of  transverse  instabilities  are  visible,  the  size  of  which  vary  substantially  across  the  width  of  the 
channel.  These  instabilities  leave  behind  highly  irregular  cellular  patterns  in  the  numerical  smoke  foil  shown 
in  Figure  3.37c.  The  initial  cellular  structures  are  relatively  small  and  regular,  but  later  complex  multi-level 
structures  appear  that  show  no  clear  distinction  between  large-scale  cells  caused  by  the  large  primary  triple 
points  on  the  reaction  front  and  small-scale  cells  caused  by  the  smaller  secondary  triple  points  that  form 
between  the  primary  instabilities.  A  close-up  view  of  the  cellular  structure  is  shown  in  Figure  3.38a.  This  is 
compared  with  a  portion  of  an  experimental  smoke  foil  of  roughly  the  same  dimensions  obtained  for  a  near- 
stoichiometric  metlrane-air  mixture  [14]  (Figure  3.38b).  The  multilevel  cellular  structures  are  also  found  on 
the  experimental  smoke  foil,  but  the  cells  in  the  photograph  seem  to  be  more  complete  than  those  found  on 
the  numerical  smoke  foil.  Similar  multilevel  instabilities  have  also  been  computed  using  a  single-component 
reaction  model  [10]  and  seem  to  be  unique  to  mixtures  with  very  large  effective  activation  energies. 

Due  to  the  irregularity  of  the  cells  produced  by  this  high-activation  energy  mixture,  there  is  no  noticeable 
difference  between  the  structures  shown  in  Figure  3.37c  for  </>  =  1  and  those  obtained  for  cf>  =  0.9  and  1.1. 
For  the  rich  ( <fi  —  1.1)  mixture,  Dqj  and  A  do  not  vary  noticeably  from  those  in  a  stoichiometric  mixture. 

3.4.3.  Detonations  in  a  Concentration  Gradient 

We  next  consider  nonunifonn  mixtures  of  methane  and  air  in  the  same  two-dimensional  channel  configuration. 
The  mixture  at  the  top  of  the  channel  ( y  =  32  cm)  is  taken  to  be  fuel  rich,  and  that  at  the  bottom  of  the 
channel  (y  =  0  cm)  is  taken  to  be  fuel  lean.  Between  these  two  extremes,  (f)  varies  linearly  from  <f)r  to 
cf>i.  Three  different  cases  are  investigated:  (1)  </>;  =  0.7  and  (f>r  =  1.3,  (2)  <f>i  =  0.8  and  </>r  =  1.2,  and  (3) 
4>i  =  0.9  and  <j>r  =  1.1.  The  behavior  of  the  detonation  is  different  for  each  of  the  three  cases,  the  reaction 
wave  propagation  velocities  of  which  are  shown  in  Figure  3.40.  For  case  1,  the  propagation  velocity  abruptly 
changes  from  ss  1700  m/s  to  «  1000  m/s  in  a  manner  similar  to  that  shown  in  Figure  3.35  for  the  </>  =  0.7 
case.  For  cases  2  and  3,  the  detonations  are  able  to  survive  and  propagate  at  speeds  between  Dcj{4n)  and 
Dcj{(t>r)-  The  same  basic  reaction  zone  structure  is  formed  for  each  of  the  three  cases  and  is  shown  in 
Figure  3.39.  The  leading  edge  of  the  reaction  zone  is  an  unstable  detonation  front.  Away  form  the  center 
of  the  channel,  where  </>  =  1,  excess  fuel  and  oxidizer  pass  through  the  propagating  detonation  wave.  These 
reactants  mix  downstream  and  burn  in  a  turbulent  diffusion  flame  that  trails  behind  the  leading  detonation. 
We  note  here  that  the  burning  in  this  diffusion  zone  is  only  qualitative.  Since  we  are  neglecting  the  diffusive 
terms  in  the  reactive  Navier-Stokes  equations,  burning  is  caused  by  numerical  diffusion  in  the  integration 
routine.  A  large  degree  of  turbulent  mixing  occurs  in  this  region,  which  also  assists  with  keeping  the  trailing 
diffusion  flame  active. 
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We  begin  by  considering  case  1,  for  which  the  variation  in  composition  d<j>  =  (f)r  —  4>i  =  0.6.  The  results 
shown  in  the  previous  section  indicate  that  a  detonation  cannot  propagate  in  a  mixture  with  <f>  =  0.7  or 
0.8  for  the  channel  size  being  used.  The  numerical  smoke  foil  shown  in  Figure  3.41a  shows  how  the  cellular 
pattern  caused  by  the  propagating  detonation  wave  in  the  case  1  mixture  evolves  from  detonation  initiation 
to  the  final  decay.  Initially,  the  detonation  is  overdriven,  allowing  a  multitude  of  smaller  cells  to  form  within 
the  channel.  Throughout  the  first  two-thirds  of  the  the  smoke  foil,  there  is  a  noticeable  difference  in  the  sizes 
of  the  cells  formed  near  the  top  of  the  channel,  where  the  mixture  is  rich,  and  those  formed  near  the  bottom, 
where  the  mixture  is  lean.  The  size  of  cells  for  the  lean  equivalence  ratios  encountered  in  case  1  are  known 
to  be  substantially  larger  than  those  for  the  rich  mixtures  encountered  in  case  1.  Hence,  the  appearance 
of  larger  cells  near  the  bottom  of  the  channel  is  physically  consistent  with  experimental  observations  in 
hydrogen-oxygen  mixtures  [89].  As  the  detonation  slows  to  a  steady  propagation  speed,  the  cell  sizes  grow 
throughout  the  channel,  and  just  prior  to  the  detonation  failure  there  is  little  difference  between  the  sizes  of 
cells  near  the  top  of  the  channel  and  those  near  the  bottom  of  the  channel.  Failure  occurs  as  triple  points 
propagate  from  the  more  reactive  mixtures  at  the  top  of  the  channel  to  the  less  reactive  mixtures  near  the 
bottom  and  all  die  before  they  can  reflect  off  the  bottom  wall. 

The  smoke  foil  for  case  2  ( <f>i  =  0.8,  </>r  =  1.2)  shows  a  system  in  which  the  survival  of  the  detonation 
is  tenuous  (Figure  3.41b).  Although  it  has  managed  to  propagate  for  a  distance  equivalent  to  10  channel 
widths,  there  is  no  guarantee  that  it  will  continue  to  propagate  indefinitely.  For  this  case  strong  triple  points 
that  originate  in  the  more  reactive  mixture  near  the  top  of  the  channel  propagate  into  the  less  reactive 
mixture  near  the  bottom  of  the  channel.  Sometimes  these  triple  points  die  out,  but  some  can  survive  and 
be  reflected  back  into  the  system.  For  some  portions  of  the  smoke  foil,  there  is  a  single  triple  point  that 
propagates  between  the  two  channel  walls.  Here,  we  have  an  added  level  of  complexity  in  that  the  reactivity 
of  the  mixture  changes  across  the  channel,  and  chance  fluctuations  in  the  propagation  speed  of  the  detonation 
could  be  sufficient  to  cause  the  triple  point  to  die.  This  behavior  causes  drastic  distortions  of  the  detonation 
surface  as  large  regions  of  reaction  zone  decouple  from  the  leading  shock  at  any  particular  instant. 

Case  3  (c f>i  =  0.9,  </>r  =  1.1)  is  more  robust  than  cases  1  and  2.  The  mixture  throughout  the  channel  is 
everywhere  sufficiently  reactive  for  detonation  propagation  to  occur.  Accordingly,  there  is  less  fluctuation 
in  the  average  propagation  velocity  (Figure  3.40),  the  magnitude  of  which  seems  to  be  slightly  less  than 
the  CJ  velocity  obtained  for  any  of  the  compositions  present  within  the  gradient.  The  cellular  patterns  left 
behind  by  this  propagating  detonation,  a  portion  of  which  is  shown  in  Figure  3.41c,  are  quite  similar  to  those 
obtained  for  the  near-stoichiometric  uniform  mixtures  (cf.  Figure  3.38a).  A  hierarchy  of  cellular  structures 
can  be  found,  the  largest  of  which  appear  to  be  incomplete.  A  sequence  of  density  plots  over  a  distance 
roughly  equal  to  one  channel  width  is  shown  in  Figure  3.42.  The  detonation  wave  appears  to  behave  in  a 
manner  similar  to  that  observed  for  a  uniform  stoichiometric  mixture. 

3.4.4.  Conclusions  and  Future  Directions 

A  single-step,  two-component  reaction  model  for  nonuniform  nrethane-air  mixtures  was  developed.  Reaction 
and  transport  coefficients  were  calibrated  to  reproduce  properties  of  detonations  and  laminar  flames  over  a 
range  of  unburned  mixture  compositions.  Chapman-Jouget  detonation  velocities  and  reaction  zone  widths 
were  found  to  match  closely  exact  values  within  this  limited  range  of  (j).  This  model  was  used  to  compute 
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the  propagation  of  detonation  waves  in  a  two-dimensional  channel  filled  with  various  uniform  methane-air 
mixtures  and  several  mixtures  in  which  the  fuel  concentration  varied  across  the  channel  width. 

The  two  leanest  uniform  mixtures  that  were  considered  (4>  =  0.7  and  0.8)  did  not  support  detonation 
propagation.  After  relaxing  to  near-CJ  velocities,  instabilities  in  the  detonation  front  caused  the  reaction 
zone  to  globally  decouple  from  the  leading  shock  wave.  Mixtures  closer  to  stoichiometric  (<j>  =  0.9,  1,  and 
1.1)  were  able  to  support  a  detonation  in  the  32  cm  channel.  Propagation  velocities  for  these  three  cases 
were  similar  and,  on  average,  were  very  close  to  their  respective  CJ  velocities.  The  trajectories  of  triple 
points  were  mapped  using  a  numerical  smoke  foil  technique.  The  cellular  structures  formed  for  these  three 
mixtures  were  nearly  identical.  Large-scale  cells  were  found  to  be  incomplete  and  composed  of  a  hierarchy 
of  smaller  cellular  structures,  which  is  characteristic  behavior  for  extremely  high  activation-energy  mixtures. 

The  reaction  zone  structure  formed  by  a  detonation  propagating  through  nonuniform  methane-air  mix¬ 
tures  in  which  the  gradient  in  fuel  concentration  is  normal  to  the  direction  of  propagation  was  also  investi¬ 
gated  for  three  different  linear  concentration  gradients  centered  about  the  stoichiometric  concentration.  The 
reaction  zone  structure,  which  involved  two  distinct  burning  regimes,  was  similar  for  all  three  cases.  The 
leading  edge  of  the  reaction  zone  is  an  unstable  detonation  wave.  Behind  the  detonation,  excess  fuel  and 
oxidizer  that  are  not  consumed  by  rich  and  lean  parts  of  the  detonation  front  mix  and  burn  in  a  turbulent 
diffusion  flame.  The  resulting  reaction-zone  structure  is  qualitatively  similar  to  the  triple  flame  detonation 
wave  (TFDW)  described  by  Calhoon  and  Sinha  [91].  Transverse  waves  were  found  to  propagate  along  the 
detonation  front  in  a  manner  similar  to  that  observed  in  uniform  mixtures.  When  the  initial  detonation 
was  overdriven,  multiple  cells  were  able  to  fit  inside  the  channel,  and  their  size  varied  from  smaller  to  larger 
between  the  regions  with  more  and  less  reactive  mixtures,  respectively,  which  is  similar  to  the  behavior 
observed  in  analogous  experiments  [89]. 

The  shape  of  the  detonation  front  was  similar  for  all  three  gradients  considered,  and  did  not  exhibit  the 
curved  structure  observed  in  experiments  [89,  90].  This  is  most  likely  due  to  the  small  system  size  used  in 
these  calculations.  The  distortions  in  the  detonation  front  caused  by  the  transverse  instabilities  were  large 
compared  to  the  width  of  the  channel,  and  local  variations  in  reaction  front  velocity  were  large  compared 
to  differences  in  CJ  velocities  over  the  range  of  (j>  considered.  Thus,  any  observed  variations  in  the  shape  of 
the  reaction  front  are  dominated  by  detonation  instabilities.  For  larger  systems  compared  to  the  detonation 
cell  size  (such  as  those  considered  in  experiments  [89,  90]),  where  many  triple  points  can  form  and  interact, 
the  influence  of  channel  size  would  be  less  pronounced.  Propagation  velocity  would  be  slower  in  the  very 
lean  and  very  rich  portions  of  the  gradient  and  would  lead  to  a  curved  detonation.  In  future  work,  we  plan 
to  perform  a  more  systematic  study  of  this  system-size  effect  as  well  as  to  explore  the  effect  of  the  various 
properties  of  the  gradient  on  detonation  propagation.  Two  important  questions  are  how  does  the  magnitude 
of  the  gradient  affect  the  detonation  cellular  structures  and  does  the  presence  of  the  turbulent  diffusion  flame 
help  to  support  detonation  propagation. 

More  investigation  is  needed  in  order  to  discuss  how  characteristics  of  the  concentration  gradient  affect 
the  structure,  propagation  speed,  and  cellular  instabilities  of  two-dimensional  detonations,  and,  likewise, 
how  the  concentration  gradient  affects  turbulent  flame  acceleration  and  DDT. 
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Figure  3.28.  Curves  of  constant  Xd  and  Si  in  the  Ea-A  parameter  space  computed  using  ID  steady-state 
flame  and  ZND  detonation  codes  in  conjunction  with  the  two-component  reaction  model  with  qM / RTq  =  39.0 
and  7  =  1.197.  The  two  curves  for  Xd  with  bold  dash-dot  lines  indicate  uncertainty  in  correlation  between 
detonation  cell  size  and  Xd ■  A  third  dash-dot  line  shows  constant  Xd  =  0.229  cm,  the  value  used  for  calibration 
of  the  single-component  reaction  model.  The  two  solid  curves  for  Si  are  computed  assuming  n0  =  6.25  x  10~6 
g/s-cm-K0  7.  The  dashed  curve  for  Si  =  38  cm/s  is  computed  for  a  higher  kd  =  3.25  x  10~5  g/s-cm-K07. 


Figure  3.29.  One-dimensional  detonation  structures  computed  using  the  single-component  reaction  model 
with  A  =  1.64  x  1013  cm3/g-s  (red  line)  and  the  two-conrponent  reaction  model  with  A  =  1.78  x  1013  cm3/g-s 
(green  line).  The  same  q,  7,  and  Ea  were  used  for  both  calculations  (cf.  Table  2  and  Table  7).  In  both  cases, 
Xd  =  0.229  cm. 
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Figure  3.30.  One-dimensional  laminar  flame  structures  computed  for  various  A,  Ea  parameter  sets.  Dashed 
lines  indicate  reaction  rates  and  correspond  to  the  right  y-axis.  Solid  lines  indicate  fuel  mass  fractions  and 
correspond  to  the  left  y-axis.  Red  lines  represent  the  flame  structure  obtained  using  a  single-component 
reaction  model  and  the  parameters  in  Table  2.  Flame  structures  computed  using  two-component  models 
with  the  same  q ,  7,  and  Ea  are  shown  for  which  Si  =  40.35  cm/s  and  Xd  =  0.229  cm  (black  lines)  and  for 
which  Si  =  40.35  cm/s  and  Xd  =  0.05  cm  (blue  lines).  For  the  former  case,  k0  was  taken  to  be  larger  than 
that  of  air  in  order  to  match  Si  and  Xd  with  the  same  set  of  parameters. 


Figure  3.31.  Velocity  of  the  leading  edge  of  the  flame  front  as  a  function  of  position  in  Configuration  17.4 
with  £  =  0.3  computed  using  the  three  reaction  models  described  in  Figure  3.30.  Experimental  data  for  two 
9.5%  methane-air  mixtures  are  taken  from  [72] 
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Figure  3.32.  Values  of  (a)  A *,  (b)  q*M/RT0:  (c)  7,  (d)  Ea/RT0,  and  (e)  kq  as  a  function  of  <f>  computed  using 
Eqs.  (3.17-3.21)  (dashed  blue  lines),  Eqs.  (3.22-3.26)  (dash-dot-dot  black  lines),  and  constant  parameters 
(dash-dot  green  lines). 
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(a)  (b) 


Figure  3.33.  (a)  Dqj  and  (b)  Xd  computed  using  one-dimensional  ZND  equations  assuming  constant  reac¬ 
tion  and  transport  parameters  (dash-dotted  green  lines),  Eqs.  (3.17-3.21)  (dashed  blue  lines),  and  Eqs.  (3.22- 
3.26)  (dash-dot-dot  black  lines)  compared  to  target  values  (solid  red  lines). 


1 024  cm 


Figure  3.34.  Schematic  of  computational  domain  for  detonation  propagation  through  a  gradient  in  fuel 
concentration.  The  equivalence  ratio  of  the  mixture  varies  linearly  from  a  maximum  of  <pr  at  the  top  of  the 
channel  to  a  minimum  of  (f>i  at  the  bottom  of  the  channel.  A  detonation  is  ignited  directly  by  a  strong  shock 
with  Mach  number  Ms  =  10  near  the  closed  end  of  the  channel.  For  all  cases,  the  width  of  the  channel  is 
32  cm,  and  its  length  is  1024  cm. 


82 


Figure  3.35.  Velocity  of  the  leading  edge  of  the  reaction  front  as  a  function  of  its  distance  from  the  closed 
end  of  the  channel  for  uniform  (0/  =  <f>r)  methane  -air  mixtures  with  <j>  =  0.7  (red  lines),  0.8  (green  lines), 
0.9  (blue  lines),  1.0  (black  lines),  and  1.1  (purple  lines)  along  with  the  value  of  Dcj  for  each  mixture  (thin 
horizontal  lines  of  the  same  color). 
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Figure  3.36.  Selected  density  maps  of  the  full  channel  width  in  the  vicinity  of  the  reaction  front  at  several 
instants  in  time  for  a  detonation  propagating  through  a  uniform  (a)  <fi  =  0.7  and  (b)  <j>  =  1.0  mixture. 
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Figure  3.37.  Numerical  smoke  foils  produced  by  detonations  propagating  through  uniform  mixtures  with 
(a)  </>  =  0.7,  (b)  </>  =  0.8,  and  (c)  (j)  =  1.0. 


(a)  (b) 


Figure  3.38.  Close-up  view  of  cellular  structures  recorded  on  (a)  a  numerical  smoke  foil  produced  by  a  det¬ 
onation  propagating  through  a  uniform  stoichiometric  (0  =  1)  mixture  and  (b)  a  portion  of  an  experimental 
smoke  foil  generated  by  detonation  of  a  uniform  near-stoichiometric  (10%  methane)  mixture.  Experiments 
were  conducted  in  a  1.05  m  diameter  and  73  m  long,  partially-closed  tube  with  smoke  foil  attached  at  the 
open  end.  The  dimensions  of  the  numerical  smoke  foil  and  the  portion  of  the  experimental  smoke  foil  shown 
in  (b)  are  roughly  the  same  32  cm  wide  by  48  cm  long.  The  direction  of  detonation  propagation  was  from 
left  to  right  in  the  figure.  Photo  courtesy  of  R.  Karl  Zipf  (NIOSH). 
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Figure  3.39.  Maps  of  the  logarithm  of  the  reaction  rate  for  a  detonation  propagating  through  a  32cm- wide 
channel  with  a  gradient  in  fuel  concentration  that  varies  linearly  from  <fir  =  1.1  to  <pi  =  0.9. 


Figure  3.40.  Velocity  of  the  leading  edge  of  the  reaction  front  as  a  function  of  its  distance  from  the  closed 
end  of  the  channel  for  methane  -air  mixtures  with  concentration  gradients  <pr  =  1.3  and  4>i  =  0.7  (red  line), 
(j)r  =  1.2  and  (j)i  =  0.8  (green  line),  and  cj)r  =  1.1  and  4>i  =  0.9  (blue  line). 
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Figure  3.41.  Numerical  smoke  foils  produced  by  detonations  propagating  through  concentration  gradients 
of  (a)  cf>r  =  1.3  and  <j>i  =  0.7,  (b)  <f>r  =  1.2  and  cj>i  =  0.8,  and  (c)  </>r  =  1.1  and  (f>i  =  0.9. 
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Figure  3.42.  Selected  density  maps  of  the  full  channel  width  in  the  vicinity  of  the  detonation  front  at 
several  instants  in  time  as  it  propagates  through  a  concentration  gradient  with  <pr  =  1.1  and  </>;  =  0.9. 


4.  Measuring  Detonability 


4.1  Background 

This  work  was  initiated  to  address  safety  concerns  related  to  explosive  gas  mixtures  that  form  in  coal  mines. 
These  are  essentially  mixtures  of  air  and  natural  gas  (NG)  released  from  coal  beds.  Similar  explosive  mixtures 
may  form  during  production,  transportation,  or  utilization  of  NG.  Though  accidental  NG  explosions  occur 
in  deflagration  regimes,  the  worst-case  scenario  involves  detonations  that  can  be  extremely  destructive  and 
generate  dynamic  pressures  up  to  10  MPa  on  reflections.  It  is  therefore  very  important  to  know  at  what 
conditions  NG-air  mixtures  can  or  cannot  detonate.  These  conditions  include  (1)  the  composition  of  the 
mixture,  (2)  the  pressure  and  temperature  of  the  mixture,  (3)  the  size  and  the  geometry  of  confinement,  and 
(4)  the  method  of  detonation  initiation. 

Here  we  will  largely  ignore  the  condition  (2)  since  we  only  consider  NG-air  mixtures  at  0.1  MPa  and 
temperatures  around  298  K,  within  seasonal  outdoor  variations.  The  remaining  three  conditions  may  vary, 
and  we  will  discuss  them  below.  The  DDT  as  a  method  for  detonation  initiation  will  be  discussed  separately 
in  section  4.1.3. 

4.1.1  Effect  of  NG  Composition  on  Detonability  of  NG-air  Mixtures 

The  composition  of  NG  is  not  strictly  defined  and  typically  includes  82-99%  of  CH4  by  volume  depending 
on  the  gas  origin  [96,  97],  as  well  as  ethane,  propane,  and  smaller  amounts  of  even  higher  hydrocarbons.  A 
few  percent  of  nitrogen  and  CO2  can  also  be  present.  A  typical  coalbed  gas  (often  called  coalbed  methane) 
usually  contains  less  than  2%  of  ethane  and  heavier  hydrocarbons,  but  sometimes  may  include  up  to  15%  of 
C02  [98,  99], 

The  variability  of  NG  composition  results  in  large  uncertainties  in  explosive  properties  of  NG-air  mix¬ 
tures.  It  has  been  shown,  for  example,  that  small  additives  of  ethane,  propane,  and  higher  alkanes  signifi¬ 
cantly  decrease  ignition  delays  in  CH4-O2  and  CH4-air  mixtures  [96,  97,  100].  For  example,  the  induction 
length  5  of  a  ZND  detonation  in  stoichiometric  methane-ethane-air  mixtures  computed  in  [100]  can  be 
approximated  as 


S/So  =  1  -  13.52  *  x  +  167.5  *x2  -  1370  *  x3  +  5846  *  x4  -  9625  *  x5  (4.1) 

where  <5o  =  2.41  cm  is  the  induction  length  for  CH4-air,  and  x  <  0.2  is  the  volume  fraction  of  ethane  in  the 
methane-ethane  mixture.  According  to  this  equation,  8%  ethane  in  the  fuel  reduces  the  induction  length  by 
a  factor  of  2. 

The  NG  used  in  GETF  experiments  described  in  Section  4.3  contained  1.5-1. 7%  of  ethane  and  only 
trace  amounts  of  higher  hydrocarbons.  According  to  Eq.  (4.1),  1.5-1. 7%  of  ethane  in  NG  decreases  the 
induction  delay  by  a  factor  of  0.8  compared  to  a  pure  CH4. 

Ignition  delays  usually  correlate  with  the  energy  required  for  the  detonation  initiation  and  the  detonation 
cell  size,  both  of  which  characterize  the  detonability  of  the  mixture.  Thus,  the  detonability  of  NG-air  varies 
with  NG  composition,  and  can  significantly  differ  from  that  of  CH4-air.  For  example,  the  strong  effect  of 


ethane  on  the  initiation  energy  has  been  shown  in  [101]  for  stoichiometric  methane-ethane-air  mixtures. 
When  the  ethane  concentration  in  fuel  decreased  from  30%  to  10%  by  volume,  the  minimum  energy  required 
for  the  detonation  initiation  in  cylindrical  geometry  increased  exponentially  from  12.5  kJ/cm  to  25  kJ/cm. 
A  rough  extrapolation  to  lower  ethane  concentrations  shows  that  the  initiation  energy  may  increase  by 
two  orders  of  magnitude  for  pure  CH4-air.  A  similar  extrapolation  of  experimental  results  [102]  obtained 
for  CH4-O2-N2  mixtures  using  the  same  experimental  setup  gives  the  initiation  energy  for  CH4-air  around 
7.5  MJ/cm.  Both  of  these  extrapolations  are  far-fetched  and  can  be  wrong  by  an  order  of  magnitude,  but 
the  results  show  that  NG  composition  can  drastically  change  the  detonability  of  NG-air  mixtures. 

4.1.2  Detonability  of  Natural  Gas- Air  Mixtures 

Because  of  the  stability  of  the  methane  molecule,  CH^-air  is  the  most  insensitive  hydrocarbon-air  mixture 
with  respect  to  detonation  initiation.  There  are  still  no  confirmed  reports  of  successful  initiation  of  self- 
supported  detonations  in  unconfined  NG-air  or  methane-air  mixtures.  An  earlier  report  [103]  that  claimed 
detonation  initiation  using  1  kg  of  TNT  inside  a  3  m  rubber  sphere  filled  with  methane-air  was  analyzed  in 
[104]  with  a  conclusion  that  the  observed  detonation  was  most  likely  decaying.  There  is  also  some  uncertainty 
in  the  exact  mixture  composition  since  the  purity  of  methane  used  in  experiments  [103]  is  not  specified. 
Report  [105]  describes  “erratic,  uneven  detonations”  initiated  in  bags  filled  with  NG-air  mixtures  using 
1001-1020  g  of  high  explosive.  These  unsteady  detonation  regimes  were  observed  only  for  NG  concentrations 
8. 6-8. 8%,  though  tested  mixture  concentrations  ranged  from  5.2%  to  12.5%.  Extrapolation  of  experimental 
results  [104]  obtained  for  various  CH4-O2-N2  mixtures  suggests  that  detonation  initiation  in  unconfined 
methane-air  would  require  about  22  kg  of  tetryl.  Large-scale  experiments  [106]  have  shown  that  37  kg  of 
Composition  B  did  not  initiate  a  detonation  in  a  stoichiometric  methane-air  (9.5%  CH4)  in  an  unconfined 
hemispherical  geometry. 

There  are  only  two  published  works  in  which  high  explosives  were  used  to  initiate  nrethane-air  detona¬ 
tions  in  tubes.  Early  experiments  [107]  carried  out  at  the  Bureau  of  Mines  Experiment  Station  at  Pittsburgh 
have  shown  that  50  g  of  a  high  explosive  (65%  ammonium  nitrate/14.6%  TNT/20.3%  sodium  chloride)  initi¬ 
ate  a  detonation-like  wave  inside  a  30.5  cm  diameter  tube  filled  with  a  methane-air  mixture  containing  9.1% 
of  CH4.  The  wave  velocity  measured  in  the  3.65  m  long  tube  varies  between  1820  and  1950  m/s,  which  is 
slightly  higher  than  the  calculated  Dcj=1801  m/s.  The  authors  conclude  that  the  observed  regime  is  not  a 
’’true”  detonation  and  it  is  partially  supported  by  the  energy  of  the  initiator. 

A  longer,  11.2  m  tube  of  the  same  30.5  cm  diameter  was  used  in  experiments  [108],  in  which  detonations 
were  initiated  using  50  g  and  70  g  of  a  slightly  more  powerful  high  explosive  (50%  ammonium  nitrate  /  50% 
TNT).  The  author  observed  the  detonation  propagation  using  multiple  optical  windows,  and  concluded  that 
steady-state  detonations  are  possible  for  mixtures  containing  from  6.3%  to  13.5%  of  methane.  The  methane 
used  in  these  experiments  was,  in  fact,  a  natural  gas  that  contained  ~  2%  of  higher  hydrocarbons. 

The  difficulties  encountered  in  initiating  self-supporting  detonations  in  methane-air  using  high  explosives 
may  create  a  perception  that  methane-air  mixtures  can  only  detonate  under  extraordinary  circumstances  not 
seen  in  the  practical  world.  A  number  of  experiments  with  different  methods  of  detonation  initiation  have 
shown  that  this  perception  is  wrong.  Detonations  in  methane-air  have  been  ignited  in  confined  conditions 
using  booster  detonations  in  more  sensitive  gas  mixtures  [109,  110]  or  dcflagration-to-detonation  transition 
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(DDT)  in  smooth  [2,  111]  or  obstructed  [72,  112-114]  channels.  Still,  available  experimental  data  on  detona¬ 
tions  in  methane-air  mixtures  are  very  limited  because  detonation  experiments  using  these  mixtures  usually 
require  large  facilities.  Ideally,  a  detonation  tube  should  be  wide  enough  to  fit  several  detonation  cells,  and 
the  tube  length  should  be  sufficient  to  observe  a  self-supporting  quasi-steady-state  detonation. 

Detonation  cell  sizes  measured  in  a  183  cm  diameter  tube  [110]  for  9.6%  methane-air  mixture  are  about 
28  cm.  Cell  sizes  measured  in  a  52  cm  diameter  tube  [72]  vary  between  13  and  29  cm  for  the  stoichiometric 
mixture,  and  can  increase  to  the  tube  width  for  lean  (8.5%  CH4)  and  rich  (12%  CH4)  mixtures.  The  wide 
range  of  measured  cell  sizes  reflects  the  fact  that  cellular  structures  in  methane-air  are  very  irregular. 

In  smaller  tubes,  metlrane-air  detonations  can  only  propagate  in  a  spinning  mode.  This  mode,  initiated 
in  a  6.35  cm  wide  square  tube  with  a  hydrogen-oxygen  booster,  has  been  used  in  [109]  to  determine  initiation 
energies  for  methane-air  mixtures  of  different  concentrations.  The  minimum  initiation  energy  of  9  MJ/m2 
has  been  obtained  for  11%  CH4.  Lower  and  upper  cletonability  limits  determined  for  methane-air  in  the 
same  experiments  are  8%  and  14.5%  CH4.  The  smallest  tube  diameter  in  which  a  methane-air  detonation 
was  ever  observed  is  3.6  cm  [114].  Detonability  limits  reported  in  [114]  for  a  5  cm  diameter  tube  are  7.5% 
and  11.5%  CH4. 

The  range  of  methane  concentrations  and  tube  diameters  discussed  are  related  to  the  ability  of  NG-air 
mixtures  and  particular  geometries  to  support  detonation  propagation,  no  matter  how  this  detonation  is 
initiated.  For  practical  systems,  however,  it  is  often  important  to  know  whether  a  detonation  can  develop 
from  a  weak  ignition  source,  that  is,  whether  a  DDT  can  occur  for  a  particular  geometry  and  mixture 
composition. 

4.1.3  Experimental  Studies  of  DDT 

There  are  several  ways  DDT  has  been  studied  experimentally.  In  the  first  type  of  experiments,  an  energetic 
gas  is  contained  in  a  long  channel  closed  at  one  end,  and  ignited  with  a  spark  near  the  closed  end.  The 
resulting,  initially  laminar  flame  accelerates,  becomes  turbulent,  and  eventually  may  produce  a  detonation. 
(See,  e.g.,  [115]).  Two  basic  scenarios  of  DDT  were  observed:  sometimes  a  detonation  appears  inside  the 
region  containing  the  turbulent  flame,  and  sometimes  it  forms  in  the  preheated,  compressed  material  between 
a  leading  shock  wave  and  the  flame.  Experiments  [116,  117]  have  shown  that  for  these  scenarios  boundary 
layers  play  a  key  role  in  the  flame  acceleration  and  turbulence  generation. 

A  second  type  of  experiment  (based  on  the  work  of  Markstein  [118])  used  a  shock- flame  interaction  to 
create  the  turbulent  flame  more  quickly  [39,  119-121].  A  flame  was  ignited  by  sparks  at  a  distance  from  a 
reflecting  wall,  and  then  a  shock  was  released  at  the  opposite  end  of  the  tube.  The  flame  interacted  first  with 
the  incident  shock  and  later  with  the  shock  reflected  from  the  end  wall.  Shock-flame  interactions  distorted 
the  flame,  making  it  turbulent  and  increasing  the  energy-release  rate  in  the  system.  These  experiments 
showed  a  number  of  specific  DDT-related  phenomena  and  well-defined  trends  on  which  simulations  could 
be  focused.  For  example,  whether  and  where  DDT  occurred  in  the  system  depended  on  the  Mach  number 
of  the  incident  shock.  For  the  lowest  Mach  numbers,  the  reflected  shock  was  amplified  by  the  turbulent 
flame,  but  there  was  no  DDT.  For  intermediate  Mach  numbers,  DDT  occurred  between  the  flame  brush  and 
the  reflected  shock.  For  the  highest  Mach  numbers,  DDT  occurred  inside  the  flame  brush.  Whenever  DDT 
occurred,  the  transition  was  through  a  very  rapid  explosion  whose  details  could  not  be  determined  from  the 
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experiment. 

A  substantial  body  of  experimental  work  has  since  shown  that  the  presence  of  obstacles  along  the  walls 
of  the  channels  results  in  much  faster  flame  acceleration  and  transition  to  a  detonation  than  if  no  obstacles 
were  present.  Channels  with  obstacles  provide  a  convenient  environment  in  which  DDT  can  be  studied  in 
a  controlled  manner.  This  technique  was  first  used  in  experiments  [122,  123],  and  it  has  been  extensively 
applied  in  combustion  research  since  then  [60-71,  124,  125].  Today  it  has  become  a  model  for  evaluation 
of  fuel  safety  using  laboratory  experiments.  Turbulent  flames  propagating  in  channels  with  obstacles  often 
accelerate  to  supersonic  velocities,  and  may  trigger  detonations  if  the  channel  is  large  enough.  The  distance  to 
DDT  in  these  experiments  depends  on  the  mixture  reactivity,  channel  width,  and  obstacle  configurations.  A 
detailed  discussion  of  experimental  results  obtained  for  different  mixtures,  including  generalized  correlations, 
can  be  found  in  the  recent  review  [126]. 

Finally,  controlled  experiments  attempting  to  observe  unconfined  DDT  [127-  129]  showed  that  a  transi¬ 
tion  to  detonation  induced  by  turbulent  flames  in  systems  without  walls  or  obstacles  is  rather  difficult.  This 
can  be  partially  attributed  to  the  lack  of  shock-reflecting  surfaces.  It  can  also  be  partially  due  to  geomet¬ 
rical  effects  of  expansion:  shocks  proceeding  from  a  deflagration  might  be  weakened,  or  turbulence  might 
be  damped  too  much  by  the  expansion,  and  so  they  become  unable  to  precondition  the  gas  for  DDT.  Wag¬ 
ner  and  coworkers  [127]  report  experiments  in  which  deflagrations  were  forced  to  DDT  by  passing  through 
screens  of  specified  mesh  sizes.  The  screens  created  turbulence  of  the  required  scale  and  intensity.  These 
experiments  suggest  that  an  unconfined  deflagration  could  make  the  transition  to  detonation  under  the  right 
conditions. 

A  related  problem  that  was  studied  experimentally  is  initiation  of  detonations  by  turbulent  jets  [128, 
130-133].  In  these  experiments,  a  jet  of  hot  combustion  products  was  injected  into  an  unburned,  cold 
mixture.  The  turbulence  created  by  the  interaction  of  this  jet  and  the  background  gas  created  a  nonuniform, 
preconditioned  region  in  which  detonation  could  occur.  These  experiments  provide  important  information  on 
the  critical  size  of  the  region  capable  of  triggering  DDT  when  the  effects  of  reflected  shocks  and  interactions 
with  walls  are  minimal. 

The  most  extreme  case  of  unconfined  DDT  is  an  initially  laminar  flame  inside  a  very  large  volume 
of  a  reactive  mixture  with  no  pre-existing  turbulence  or  shocks.  In  the  absence  of  effective  mechanisms 
for  turbulence  generation,  such  as  shock-flame  interactions  and  wakes,  laminar  flames  become  turbulent 
through  other,  relatively  slow  mechanisms,  such  as  Darrieus-Landau  and  thermal- diffusive  instabilities.  In 
the  presence  of  gravity,  the  buoyancy-driven  Rayleigh-Taylor  (RT)  instability  might  become  a  dominant 
mechanism  for  turbulence  generation.  The  developing  turbulence,  however,  might  be  damped  by  expansion. 
The  size  of  an  unconfined  system  in  which  DDT  could  be  observed  is  expected  to  be  relatively  large.  It  has 
been  suggested,  for  example,  that  DDT  might  occur  in  very  large  vapor  clouds  [134,  135]. 

Analysis  of  experimental  data  shows  that  turbulence  plays  an  important  role  in  DDT.  Several  apparently 
different  mechanisms  for  the  DDT  have  been  described,  each  including  the  effects  of  turbulence  and  formation 
of  shocks.  On  large  scales,  turbulence  deforms  the  flame  front  and  increases  its  surface  area.  On  small  scales, 
it  can  broaden  the  flame  front  and  cause  mixing.  The  result  is  an  extended  turbulent  “flame  brush”  in  (or 
near)  which  a  series  of  explosions  occurs.  One  of  these  explosions  finally  leads  to  a  detonation.  Other  routes 
to  detonation  may  include  an  explosion  in  the  boundary  layer,  or  an  explosion  in  the  region  between  the 
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leading  shock  and  flame  brush. 

This  summary  of  observational  evidence  gives  a  general  picture  of  when  and  where  DDT  might  occur,  and 
the  importance  of  fluid  instabilities  in  creating  the  turbulent  background  in  which  a  detonation  is  more  likely 
to  appear.  Available  experimental  data,  however,  do  not  directly  show  how  a  detonation  appears.  Studies 
of  intrinsic  mechanisms  by  which  detonations  form  require  a  combination  of  experimental,  numerical,  and 
theoretical  approaches  [5]. 

4.2  Gas  Explosion  Test  Facility  at  NIOSH  Lake  Lynn  Laboratory 

Experiments  were  carried  out  at  the  Gas  Explosion  Test  Facility  (GETF)  at  NIOSH  Lake  Lynn  Laboratory 
(LLL).  The  facility  consists  of  the  detonation  tube,  supporting  systems  for  gas  injection,  mixing,  sampling 
and  analysis,  the  ignition  system,  the  diagnostic  equipment,  and  a  remote  control  system.  Pictures  of  the 
facility  are  shown  in  Figure  4.1,  and  a  schematic  diagram  is  shown  in  Figure  4.2. 

4.2.1  Detonation  Tube 

The  tube  is  located  in  the  upper  quarry  of  LLL,  and  is  fabricated  from  sections  of  hot-rolled  and  welded 
industrial  steel  pipe.  The  total  length  of  the  tube  is  73  m,  and  the  internal  diameter  is  105  cm.  The  straight 
tube  rests  in  plane  on  18  equally-spaced  concrete  blocks  weighing  about  1.8  t  each,  and  is  attached  to  these 
blocks  with  9.5  mm  galvanized  steel  cable.  The  tube  is  inclined  about  2  deg,  so  that  the  open  end  is  about 
2.5  m  higher  than  the  closed  end. 

The  closed  end  is  permanently  sealed  with  a  12.7  mm  welded  steel  plate,  and  rests  against  a  14  t  concrete 
block  poured  between  the  tube  and  a  5-m-high  earthen  berm,  as  shown  in  Figure  4.1  (left).  The  total  weight 
of  the  pipe,  concrete  supports,  and  the  earthen  berm  that  resists  the  tube  recoil  is  at  least  175  t.  The 
open  end  of  the  tube  is  located  about  15  m  from  the  20-m-high  quarry  wall.  A  15  t  concave  ramp  deflector 
fabricated  from  curved  steel  H-beams  is  installed  at  about  2  m  from  the  open  end,  leaning  against  another 
berm.  The  deflector  is  designed  to  direct  shock  waves  upward  into  the  atmosphere  in  order  to  minimize 
strong  shock  reflections  from  the  quarry  walls. 

The  steel  pipe  has  a  wall  thickness  of  9.5  mm,  minimum  yield  strength  of  248  MPa,  and  an  ultimate 
tensile  strength  of  414  MPa.  For  stoichiometric  methane-air  mixtures,  the  quasi-static  Chapman-Jouguet 
detonation  wave  pressure  is  1.66  MPa,  and  that  pressure  could  persist  for  time  periods  longer  than  0.1  second. 
Based  on  hoop  stress  calculations  with  this  quasi-static  pressure,  the  safety  factor  against  deformation  is 
about  2.70,  and  the  safety  factor  against  rupture  is  4.50.  Higher  pressure  pulses  with  shorter  durations  affect 
the  tube  walls  in  the  vicinity  of  the  detonation  front,  during  shock  reflections,  and  during  the  detonation 
initiation.  The  highest  dynamic  pressure  recorded  so  far  at  GETF  is  about  7  MPa,  and  the  tube  shows  no 
evidence  of  deformation  or  damage. 
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4.2.2.  Gas  Mixing  and  Analysis 

A  very  important  part  of  GETF  is  the  gas  mixing  system  that  enables  creation  of  homogeneous  methane-air 
test  mixtures  inside  the  tube  with  predefined  equivalence  ratios.  The  mixing  system  consists  of  an  enclosed 
blower  fan,  an  intake  manifold,  and  an  exhaust  manifold.  The  2.25  kw  blower  fan  delivers  26.5  m3/min 
operating  at  27.2  cm  of  water  static  pressure.  Manifolds  are  made  with  10  cm  diameter  pipe,  and  they  are 
connected  to  the  detonation  tube  through  7.5  cm  inlet  and  outlet  riser  pipes.  The  riser  pipes  are  equipped 
with  manual  10  cm  ball  valves  that  can  be  used  to  balance  the  flow  from  each  inlet  and  outlet;  however,  in 
these  experiments,  all  the  valves  were  open.  There  are  also  two  remotely  operated  ball  valves  between  the 
blower  fan  and  manifolds.  These  are  open  during  the  mixing  process,  but  closed  before  ignition  to  protect 
the  fan. 

Prior  to  filling  the  explosion  tube  with  methane  to  obtain  the  desired  test  mixture  composition,  a  0.15- 
mm-thick  plastic  diaphragm  is  placed  over  the  open  end  of  the  tube.  Natural  gas  is  introduced  into  the 
intake  manifold  at  the  beginning  of  the  mixing  process. 

As  shown  in  the  schematic  in  Figure  4.2,  the  intake  manifold  has  six  evenly-spaced  inlets,  with  three 
located  on  each  side  of  the  blower  fan.  The  exhaust  manifold  has  seven  outlets,  with  three  on  one  side  of  the 
blower  fan  and  four  on  the  other.  The  first  riser  pipe  located  about  100  cm  from  the  closed  end  of  the  tube 
is  an  outlet,  that  is  mixture  flows  out  of  the  manifold  into  the  tube.  The  second  riser  is  an  inlet,  in  which 
mixture  flows  into  the  manifold  from  the  tube.  The  third  riser  is  an  outlet;  the  fourth  riser  is  an  inlet,  and 
so  on.  The  last  riser  located  about  82  cm  from  the  open  end,  is  an  outlet. 

When  the  blower  fan  operates,  the  gas  mixture  in  the  tube  is  drawn  into  the  inlet  risers,  through  the 
intake  manifold,  through  the  blower,  through  the  exhaust  manifold,  and  then  reinjected  through  the  outlet 
risers.  Since  both  inlets  and  outlets  are  distributed  along  the  length  of  the  detonation  tube,  the  system 
has  multiple  recirculation  circuits  that  merge  in  the  intake  manifold  and  split  in  the  exhaust  manifold,  thus 
providing  efficient  gas  mixing.  The  volume  of  the  detonation  tube  is  63.4  m3.  The  2.25  kw  blower  fan  moves 
26.5  m3/min,  and  thus  the  air  change  time  is  less  than  3  min.  To  achieve  a  homogeneous  test  mixture, 
typical  mixing  times  range  from  30  to  45  min. 

A  sample  draw  system  enables  the  remote  collection  and  analysis  of  samples  of  the  test  mixture.  During 
gas  mixing,  samples  are  drawn  continuously  from  the  closed  end,  middle,  and  open  end  of  the  tube  back 
to  an  infrared  gas  analyzer  (Siemens  Ultramat  23)  located  about  300  m  away  from  the  detonation  tube 
in  the  control  building.  Collected  samples  are  then  returned  to  the  detonation  tube  by  a  reinjection  line, 
except  for  a  small  portion  passing  through  the  gas  analyzer.  The  sample  draw  system  provides  real-time  gas 
concentrations  at  three  different  locations  along  the  tube  and  assures  homogeneity  of  the  test  mixture. 

The  mixture  composition  in  the  tube  is  controlled  in  three  different  ways.  First,  a  pre-calculated  amount 
of  natural  gas  is  injected  into  the  tube  using  a  conventional  natural  gas  meter  (Sensus  135  S-275).  Second, 
gas  samples  are  analyzed  with  the  remote  infrared  gas  analyzer  during  mixing.  Finally,  two  samples  of  the 
gas  mixture  are  drawn  after  the  mixing  is  completed,  and  later  analyzed  off-site  using  gas  chromatography. 
In  most  experiments  with  this  detonation  tube,  the  desired  concentration  of  natural  gas  in  the  test  mixture 
was  created  with  the  accuracy  ±0.2%  by  volume. 
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4.2.3.  Ignition 

To  ignite  a  detonation  in  the  NG-air  test  mixture,  a  large  plastic  bag  filled  with  a  stoichiometric  methane- 
oxygen  mixture  is  placed  at  the  closed  end  of  the  tube.  The  bag  is  made  of  0.15  mm  plastic  film  and  has  a 
cylindrical  shape  with  a  diameter  of  about  about  1  m.  The  bag  volume  typically  varies  between  about  3  m3 
and  6  nr3  from  one  series  of  experiments  to  another.  Smaller  bags,  1.5,  0.75,  and  0.375  m3,  were  also  used 
in  several  experiments  to  test  the  effect  of  the  ignition  energy. 

To  inflate  the  bag,  separate  methane  and  oxygen  tanks  are  used  to  meter  pre-calculated  amounts  of 
each  gas  via  orifices  into  the  bag.  Gas  was  administered  through  separate  lines  with  check  valves  and  flash 
back  arrestors  in  each  line.  The  bag  is  inflated  remotely  after  NG-air  mixing  in  the  tube  is  complete. 

Detonation  in  the  stoichiometric  methane-oxygen  mixture  is  easily  initiated  with  a  non-electric  ^8 
blasting  cap  containing  about  0.45  g  pentaerythritol  tetranitrate  (PETN)  and  placed  at  about  0.5  m  from 
the  closed  end  of  the  tube.  The  methane-oxygen  mixture  has  a  detonation  cell  size  of  2-3  mm,  propagates 
with  a  velocity  close  to  D cj  =  2490  m/s,  burns  all  the  mixture  in  the  bag,  and  generates  a  strong  shock  in 
the  tube.  This  shock  may  or  may  not  ignite  a  detonation  in  the  NG-air  test  mixture,  depending  on  the  test 
mixture  composition  and  the  size  of  the  igniter  bag. 

In  experiments  without  the  methane-oxygen  booster  bag,  the  NG-air  test  mixture  was  ignited  directly 
by  an  electric  match  placed  near  the  closed  end  of  the  tube.  This  type  of  ignition  produced  a  flame. 

4.2.4.  Baffles 

For  DDT  experiments,  for  which  there  was  no  booster  bag,  orifice  plates  (or  baffles)  with  blockage  ratio  of 
BR=0.6  were  placed  inside  the  tube  1  m  apart.  The  baffles  are  made  of  12.8  mm  thick  steel  and  connected 
with  steel  rods  as  shown  in  Figure  4.3. 

4.2.5.  Diagnostics 

The  tube  is  equipped  with  23  quartz-piezoelectric-type  pressure  transducers  and  23  light  sensors  placed  in 
pairs  (pressure  +  light)  approximately  every  3  m  on  the  tube  wall,  as  shown  in  Figure  4.2.  The  sensors 
are  connected  to  two  separate  24-channel  data  acquisition  units,  one  of  which  records  all  pressure  signals, 
another  records  all  light  signals,  and  both  units  also  record  a  synchronization  signal  from  a  signal  generator. 

The  light  sensors  are  silicon  phototransistors  (Optek  NPN  OP800WSL)  with  a  rise  time  of  about  2  gs 
and  a  spectral  response  from  600  to  1100  nrn.  Sensors  are  mounted  at  the  end  of  a  small  tube  that  is  15  cm 
long  with  an  inside  diameter  of  0.32  cm.  The  Held  of  view  for  the  light  sensor  is  about  5  cm  on  the  opposite 
side  of  the  105-cm-diameter  detonation  tube. 

The  piezoelectric  pressure  sensors  (PCB  Piezotronics  CA102B04)  have  a  range  of  0-6900  kPa  and  1  gs 
rise  time.  In  early  experiments  with  the  detonation  tube,  the  pressure  sensors  were  threaded  directly  into 
the  detonation  tube  steel.  The  speed  of  sound  in  mild  steel  is  about  5900  m/s,  and  extraneous  vibrations 
induced  in  the  steel  by  the  methane-air  combustion  arrive  at  the  pressure  sensor  prior  to  the  arrival  of  the 
detonation  wave,  which  travels  on  the  order  of  1800  m/s.  To  reduce  this  noise,  pressure  sensors  are  decoupled 
from  the  tube  using  double-studded  silicone  gel  anti-vibration  mounts,  as  shown  in  Figure  4.4. 

The  data  acquisition  system  (DAS)  consists  of  two  separate  24-clrannel  units  (National  Instruments 
CompactRIO)  with  a  sampling  frequency  of  50  kHz.  This  frequency  limits  the  time  resolution  for  all  recorded 
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signals  to  20  /is.  The  time  error  between  two  signals  generated  by  a  wave  propagating  between  two  sensors 
located  3  m  apart  can  reach  40  frs.  For  the  detonation  velocity  1800  m/s,  the  uncertainty  in  position  of  the 
detonation  wave  is  0.072  m  in  3  m  or  2.4%.  Therefore,  the  uncertainty  in  the  measured  detonation  velocity 
is  also  2.4%  or  ±43  m/s.  If  the  velocity  is  constant,  the  uncertainty  is  significantly  reduced  using  multiple 
sensors  at  greater  distances  apart. 

Before  each  experiment,  two  smoke  foils  are  placed  between  0.15  and  2.8  m  from  the  open  end  of  the  tube. 
The  1.2-rn  by  1.2-m  foils  are  made  of  1.6-mm-thick  aluminum  sheets  rolled  to  a  0.5-rn  radius  of  curvature, 
and  secured  to  the  inside  top  of  the  tube  with  bolts  around  the  perimeter  of  the  foils.  The  upstream  edge 
of  each  foil  is  inserted  into  a  welded  metal  slot.  The  foils  are  covered  with  soot  produced  by  acetone  flames. 

4.2.6.  Remote  Control  System 

Starting  from  the  injection  of  NG  into  the  detonation  tube,  all  operations  necessary  for  the  preparation  and 
conduction  of  experiments  are  controlled  remotely  from  the  LLL  control  building  located  about  300  m  from 
the  tube.  Firing  the  primary  ignition  circuit  triggers  the  data  acquisition  system.  In  case  the  primary  circuit 
fails  to  ignite  the  methane-oxygen  bag  and  the  test  mixture,  two  additional  independent  ignition  circuits  are 
available  to  ignite  the  flammable  mixture  using  electric  matches.  It  is  also  possible  to  vent  the  test  mixture 
to  atmosphere  using  the  gas  mixing  system. 
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Figure  4.1  Gas  Explosion  Test  Facility  (GETF)  at  NIOSH  Lake  Lynn  Laboratory.  Views  from  the  open 


(a)  and  closed  (b)  ends. 
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Figure  4.2  GETF  schematics. 
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Figure  4.3  Steel  baffles  with  blockage  ratio  BR=0.6  installed  inside  GETF  tube  for  DDT  experiments. 


Figure  4.4  Light  (a)  and  pressure  (b)  sensors  and  anti- vibration  mount  (c). 
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4.3  Interpreting  GETF  results 

GETF  was  used  for  two  types  of  experiments: 

1.  Direct  initiation  of  detonations  in  NG-air  mixtures  using  a  booster  detonation  in  methane-oxygen. 
The  results  include  detonability  limits,  detonation  cell  patterns  for  the  full  range  of  methane  con¬ 
centrations,  and  detonation  velocities  and  pressures. 

2.  Flame  acceleration  and  DDT  in  NG-air  mixtures  using  baffles.  The  mixture  was  ignited  by  an 
electrical  match.  The  results  include  flame  and  shock  velocities,  pressures,  DDT  distances,  and 
smoke  foil  records  obtained  for  various  numbers  of  baffles  and  mixture  compositions. 

Here  we  discuss  these  results  in  more  detail. 

4.3.1.  Detonations  in  NG-Air  Mixtures 

During  2009  and  2010,  more  than  40  explosion  tests  were  conducted  to  determine  the  detonation  charac¬ 
teristics  of  various  mixtures  ranging  from  4  to  19%  NG  in  air.  Each  experiment  attempted  direct  initiation 
of  detonation  in  the  test  mixture  using  a  bag  containing  about  2.9  m3  of  methane-oxygen  mixture  as  an 
initiator.  The  diagnostic  equipment  used  in  these  experiments  allowed  for  the  identification  of  detonations 
using  four  methods:  1)  comparison  of  measured  flame  and  shock  velocities  to  the  computed  D cj  values 
derived  from  data  taken  by  pressure  and  light  sensors  placed  in  pairs  every  3  m  on  the  tube  wall;  2)  com¬ 
parison  of  the  recorded  pressures  to  the  computed  pressure  profiles  for  ZND  detonations;  3)  measurement  of 
the  separation  between  the  shock  and  the  flame  using  the  time  difference  between  pressure  and  light  signals 
recorded  at  the  same  location;  and  4)  the  patterns  produced  on  smoke  foils  placed  near  the  open  end  of  the 
tube. 

The  results  of  direct  initiation  experiments  are  summarized  in  Figure  4.5,  which  shows  the  measured 
flame  velocity  versus  methane  concentration.  The  theoretical  CJ  detonation  velocity  computed  with  CHEE¬ 
TAH  is  also  included  in  Figure  4.5  as  a  solid  line.  Detonation  in  the  test  mixture  is  sustainable  for  methane 
concentrations  from  about  5.3  to  about  15.5%  by  volume.  The  difference  between  the  measured  and  theo¬ 
retical  detonation  velocities  is  less  than  10  m/s  or  0.6%.  Beyond  the  detonation  limits,  the  shock  velocity 
decreased  to  about  500  to  800  m/s  for  mixtures  less  than  5%  methane  in  air  and  to  about  700  to  800  m/s  for 
mixtures  greater  than  16%  methane  in  air.  These  observed  shock  velocities  are  about  1/2  D cj-  Test  mix¬ 
tures  with  less  than  4.9%  or  greater  than  16.0%  of  methane  did  not  sustain  combustion,  either  by  detonation 
or  deflagration,  when  ignited  by  the  strong  methane-oxygen  ignition  source  used  in  these  experiments. 

Figure  4.6  shows  the  measured  detonation  pressure  versus  methane  concentration.  The  theoretical  CJ 
detonation  pressure  is  also  included  as  a  solid  line.  Again,  there  is  good  agreement  between  measured  and 
theoretical  values,  considering  that  the  estimates  of  the  detonation  pressure  behind  the  main  shock  are  ±0.3 
MPa.  The  maximum  CJ  pressure,  1.77  MPa,  occurs  for  methane  concentration  about  10%,  and  declines  as 
the  lower  and  upper  limits  are  approached.  Even  at  the  limits  of  5.3  and  15.5%,  CJ  pressures  are  still  high 
at  about  1.2  MPa  and  1.5  MPa,  respectively.  Thus,  a  methane-air  mixture  develops  variable  but  substantial 
pressures  over  the  entire  detonable  range. 

Figure  4.7  presents  measured  detonation  cell  sizes  as  a  function  of  mixture  composition.  Tests  for 
mixture  compositions  in  the  range  of  5.3  to  15.5%  CH4  produced  cell  patterns,  and  those  outside  the  range 
did  not.  The  number  of  cells  in  a  foil  ranged  from  1  to  8,  depending  on  the  mixture.  Near-stoichiometric 
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mixtures  produced  a  detonation  cell  size  of  about  25  cm.  With  a  7%  test  mixture,  cell  size  increased  to 
about  55  cm,  and  with  a  15.4%  mixture,  a  single  cell  measuring  105-cm-wide  by  170-cm-long  was  recorded. 
The  variability  of  the  cell  measurements  is  on  the  order  of  30%  due  to  the  irregularity  of  the  cell  patterns. 
Figure  4.7  also  includes  the  computed  cell  size  based  on  a  model  by  Gavrikov  et  al.  [59].  There  is  a  good 
agreement  between  the  measured  and  theoretical  cell  sizes  considering  the  high  variability  of  the  cell  size 
measurements . 

Cell  patterns  left  on  smoke  foils  by  developed  detonations  were  very  irregular  and  showed  secondary  cell 
structures  inside  primary  cells.  Figures  4.8-4.10  show  examples  of  experimental  and  computed  cell  patterns. 
There  is  a  good  agreement  between  the  experimental  and  numerical  results  including  the  cell  irregularity 
(Figure  4.8),  secondary  patterns  inside  main  cells  (Figure  4.9),  and  fine  cell  structures  produced  by  transverse 
detonations  (Figure  4.10). 

The  most  important  practical  result  of  these  experiments  is  the  new  extended  range  of  detonability  for 
methane-air  mixtures  observed  in  the  large,  105  cm  diameter  tube.  The  new  range,  5.3  to  15.5%  CH4,  almost 
encompasses  the  entire  normal  combustion  limits  of  5  to  16%  methane  in  air  as  reported  by  Cashdollar  et  al. 
[136].  Figure  4.11  summarizes  detonability  limits  measured  for  nrethane-air  mixtures  by  different  research 
groups  in  tubes  of  different  diameters  d,  and  using  different  detonation  initiation  methods.  Matsui  [114] 
used  DDT  with  a  Schelkin  spiral,  Wolanski  et  al.  [109]  used  hydrogen-oxygen  booster,  Kogarko  [108]  used  a 
high  explosive,  and  Zipf  et  al.  [14]  used  methane-oxygen  booster.  Most  of  the  data  fit  nicely  on  a  bell-shaped 
curve  which  shows  wider  limits  for  larger  systems.  Similar  curves  were  obtained  by  Matsui  [114]  for  other 
hydrocarbons  and  hydrogen. 

Kuznetsov  et  al.  [113]  studied  the  flame  acceleration  in  obstructed  channels.  Though  detonations  were 
not  observed  in  these  experiments,  the  boundary  between  subsonic  and  supersonic  flames  shown  by  cross 
points  in  Figure  4.11  is  very  close  to  the  curve  that  represents  detonability  limits.  This  may  imply  that 
supersonic  flames  in  obstructed  channels  partially  rely  on  shock  compression  for  an  accelerated  energy 
release.  For  mixtures  that  cannot  detonate,  the  energy  release  does  not  accelerate  enough  in  response  to  the 
shock  compression,  and  therefore  supersonic  flames  do  not  develop. 

4.3.2.  DDT  in  NG-Air  Mixtures 

The  initial  series  of  DDT  experiments  in  GETF  show  that  detonations  in  105  cm  diameter  tube  can  develop 
from  a  weak  ignition  source.  When  9  or  more  baffles  with  blockage  ratio  BR=0.6  were  installed  at  the 
beginning  of  the  tube,  detonations  developed  in  the  smooth  section  of  the  tube  starting  at  about  50  m  from 
the  closed  end.  These  detonations  were  identified  using  pressure  and  light  sensors,  and  smoke  foils  the  same 
way  as  in  the  direct  initiation  experiments.  The  mechanism  of  DDT  in  the  smooth  section  is  different  from 
that  in  our  simulations  performed  for  the  obstructed  tube.  The  simulations  predict  that  DDT  will  occur 
in  the  obstructed  section  at  about  20  m,  that  is,  with  at  least  20  baffles  installed.  So  far,  the  experiments 
were  carried  out  with  up  to  18  baffles,  and  produced  detonations  only  in  the  smooth  section,  far  from  the 
obstructed  area.  The  pressures  at  the  end  of  the  18-baffle  section  were  very  high,  indicating  that  several 
more  baffles  may  initiate  detonations  in  the  obstructed  section.  These  pressures  lead  to  significant  damage 
of  the  last  several  baffles.  We  expect  that  a  new  series  of  experiments  with  reinforced  baffle  mounts  will 
show  the  DDT  in  the  obstructed  area. 
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Figure  4.5  Flame  velocity  measured  in  direct  initiation  experiments  (points)  as  a  function  of  mixture 
composition.  Solid  line  shows  computed  CJ  velocity. 
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Figure  4.6  Pressure  measured  behind  the  main  peak  in  direct  initiation  experiments  (points)  as  a  function 
of  mixture  composition.  Solid  line  shows  computed  CJ  pressure. 
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Figure  4.7  Detonation  cell  size  measured  in  direct  initiation  experiments,  (low,  high,  and  average  data 
points).  Solid  line  shows  computed  cell  size  based  on  a  model  by  Gavrikov  et  al.  [59]. 


Figure  4.8  Experimental  (left)  and  computed  (right)  detonation  cells.  The  fragment  of  experimental  smoke 
foil  is  25  cm  wide  and  obtained  for  10.4%  CH4. 
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Figure  4.9  Experimental  (left)  and  computed  (right)  secondary  detonation  cells.  The  fragment  of  experi¬ 
mental  smoke  foil  is  28  cm  wide  and  obtained  for  6.2%  CH4. 


Figure  4.10  Experimental  (left)  and  computed  (right)  fine  cellular  strictures  produced  by  transverse  deto¬ 
nations.  The  fragment  of  experimental  smoke  foil  is  18  cm  wide  and  obtained  for  15.5%  CH4. 
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Figure  4.11  Detonability  limits  measured  for  methane-air  mixtures  by  different  research  groups  in  tubes 
of  different  diameters  d.  Solid  line  is  the  best  fit  for  most  of  experimental  points.  Cross  points  show  the 
boundary  between  subsonic  and  supersonic  flames  in  obstructed  channels. 
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5.  Application  to  Coal  Mining 

Accidental  gas  explosions  in  underground  coal  mines  generate  pressures  that  can  result  in  the  loss  of  human 
life.  To  prevent  these  explosions,  the  operational  area  of  mines  are  constantly  ventilated  to  keep  methane 
concentrations  below  the  lower  flammability  limit.  Unventilated,  abandoned  areas  are  usually  sealed  off  with 
concrete  walls  (seals).  The  seals  are  designed  to  withstand  explosion  pressures  if  an  accidental  explosion 
occurs  in  a  sealed  area.  It  is,  therefore,  important  to  know  the  conditions  under  which  explosions  develop 
and  the  maximum  pressures  these  explosions  can  generate. 

The  main  objective  of  this  work  is  to  predict  these  maximum  pressures.  In  particular,  a  detonation  is 
considered  as  the  worst-case  explosion  scenario,  one  that  can  generate  extremely  high  pressures  ranging  from 
1.2-1. 7  MPa  at  the  Chapman-Jouguet  point  (Figure  4.6)  to  10  MPa  at  reflections.  The  ability  of  a  mixture 
of  a  reactive  gas  and  air  to  develop  a  detonation  depends  on  several  conditions  that  can  be  separated  in  two 
main  groups. 

The  first  group  defines  absolute  detonability ,  which  is  the  ability  of  a  mixture  to  support  a  detonation 
in  principle,  no  matter  how  this  detonation  is  ignited.  The  conditions  determining  absolute  detonability 
include  the  mixture  composition,  the  background  pressure  and  temperature,  and  the  system  size. 

The  second  group  defines  conditional  detonability,  which  is  the  ability  of  a  mixture  to  develop  a  deto¬ 
nation  for  a  particular  geometry  and  ignition  source.  For  weak  ignition  sources  that  do  not  initially  involve 
strong  shocks  and  detonations,  conditional  detonability  is  defined  by  the  ability  of  a  system  to  undergo  a 
deflagration-to-detonation  transition  (DDT). 

In  the  Introduction  to  this  report,  we  asked  an  overriding  question  that  is  related  both  to  the  absolute 
detonability  of  the  gas  mixture  and  conditional  detonability  of  the  particular  system: 

Given  a  large  enough  volume  of  a  flammable  mixture  of  NG  and  air,  such  as  may  exist  in  a  coal 

mine,  can  a  weak  spark  ignition  develop  into  a  detonation? 

First,  we  answer  this  overriding  question  in  part:  We  have  shown  both  computationally  and  experimentally 
that  for  systems  at  the  coal  mine  scale,  where  tunnels  measure  up  to  3  m  high,  6  m  wide  and  hundreds  of 
meters  long,  detonation  of  methane-air  mixtures  is  supportable  across  a  range  of  compositions  from  5.3  to 
15.5  %  methane  in  air.  The  least  dimension  of  coal  mine  tunnels,  the  height,  is  usually  greater  than  1  m,  and 
this  dimension  is  much  larger  than  the  detonation  cell  size,  which  is  about  20  cm  for  a  stoichiometric  mixture 
of  methane  in  air.  While  it  may  require  more  energy  and  be  more  difficult  to  develop  or  initiate  a  detonation 
in  methane-air  mixtures  near  the  detonation  limits,  it  is  possible  to  achieve  this.  Such  detonations  propagate 
in  a  stable  manner  as  long  as  a  detonable  methane-air  mixture  is  available. 

Below  we  summarize  the  results  of  this  study  in  terms  of  the  absolute  and  conditional  detonability, 
discuss  DDT,  and  return  with  a  more  general  answer  to  the  overriding  question. 

5.1.  Absolute  Detonability 

Direct  initiation  experiments  at  GETF  have  shown  that  NG-air  mixtures  containing  between  5.3  and  15.5% 
of  NG  can  detonate  in  a  tube  with  a  diameter  of  105  cm.  These  detonability  limits  become  narrower  for 
smaller  systems  (Figure  4.11),  and  are  likely  to  widen  for  larger  systems.  In  particular,  larger  detonability 
limits  are  expected  in  coal  mining  tunnels,  which  can  be  up  to  3  m  high  and  6  m  wide. 
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As  the  system  size  increases,  detonations  are  able  to  propagate  in  leaner  and  richer  mixtures.  Though 
the  energy  required  for  detonation  initiation  increases  sharply  for  mixtures  far  from  stoichiometric,  the 
chemical  energy  available  for  the  detonation  initiation  through  DDT  also  increases  with  the  system  size. 
Even  if  the  system  geometry  does  not  allow  DDT  near  the  detonation  limits,  systems  with  concentration 
gradients  may  develop  detonations  through  DDT  at  some  locations  with  near-stoichiometric  concentrations. 
These  detonations  may  then  spread  to  lean  or  rich  mixtures  in  other  locations. 

The  importance  of  the  observed  wider  detonability  limits  for  larger  systems  goes  well  beyond  large- 
scale  explosions  and  coal  mine  applications.  When  detonability  is  concerned,  the  effect  of  a  larger  scale  is 
often  similar  to  the  effect  of  a  higher  pressure,  as  recently  shown  in  experiments  [113].  In  other  words,  we 
should  expect  wider  detonability  limits  in  small  systems  at  elevated  pressures.  This  means  that  the  wider 
detonability  limits  should  be  applied  to  the  design  of  engines  that  burn  nrethane-air  mixture  at  increased 
pressures. 

5.2.  Deflagration-to-Detonation  Transition 

Even  if  the  mixture  composition  is  within  the  detonability  limits,  it  does  not  mean  that  a  detonation  will 
develop  for  a  given  geometry.  Unconfined  NG-air  clouds,  which  might  form  when  NG  leaks  from  a  vessel 
into  the  open  air,  are  very  unlikely  to  develop  detonations  from  a  weak  ignition  source.  In  this  case,  there  is 
nothing  to  contain  the  pressure  generated  by  the  propagating  flame  or  to  generate  strong  turbulence.  NG-air 
clouds  in  a  coal  mine  are  generally  not  unconfined,  except  perhaps  for  the  smallest  clouds  that  might  form  at 
an  active  face.  In  confined  spaces,  such  as  coal  mine  tunnels,  detonations  can  develop  more  easily,  especially 
when  obstructions  are  present.  DDT  in  NG-air  mixtures  was  observed  in  smooth  and  obstructed  channels 
of  various  diameters  that  range  from  3.6  cm  [114]  to  105  cm  (GETF  experiments). 

The  important  parameter  characterizing  the  ability  of  a  given  mixture  to  develop  a  detonation  in  a 
given  geometry  is  TddT)  the  run-up  distance  to  DDT  from  the  ignition  source.  This  parameter  is  closely 
related  to  the  minimum  volume  of  reactive  mixture  required  for  the  detonation  development. 

In  smooth  tubes,  flames  develop  relatively  slowly,  especially  when  ignited  by  a  weak  spark  without 
a  pre-existing  background  turbulence.  For  nrethane-air  mixtures,  this  type  of  ignition  in  smooth  tubes 
produces  a  continuous  flame  acceleration,  but  DDT  was  not  observed  when  the  tube  length  was  30  nr  and 
the  diameter  was  40  cm  [2],  and  in  GETF  experiments  without  obstacles  for  which  the  tube  length  is  73  nr 
and  the  diameter  is  105  cm.  The  observed  continuous  flame  acceleration  indicates  that  detonations  would 
most  likely  appear  in  longer  tubes,  which  means  that  L  o dt  is  greater  than  30  m,  for  d  =  40  cnr,  and  greater 
than  73  m  for  d  =  105  cm,  but  most  likely  finite. 

The  idealized  case  of  a  weak  ignition  in  a  straight  smooth  channel  without  turbulence  is  unlikely  to 
occur  in  coal  mines.  Typical  mine  tunnels  have  very  rough  walls  with  crosscut  openings  to  other  tunnels, 
can  be  partially  obstructed  by  rocks  or  mining  equipment,  and  may  have  a  pre-existing  turbulence  due  to 
the  forced  ventilation.  In  these  conditions,  the  flanre  developments  would  be  significantly  accelerated.  Even 
in  smooth  channels,  the  flame  evolution  is  accelerated  if  the  ignition  is  more  violent.  This  does  not  change 
the  mechanism  of  flame  acceleration  and  DDT,  but  accelerates  the  process  by  introducing  initial  turbulence. 

For  example,  DDT  in  smooth  tubes,  or  smooth  sections  of  partially  obstructed  tubes  was  observed  for 
a  more  violent  flame  ignition.  In  experiments  [2] ,  the  gas  was  ignited  by  a  turbulent  flame  jet  generated  by 
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combustion  in  a  separate  vessel  connected  to  the  main  tube.  In  experiments  [112],  a  high-speed  flame  jet 
was  created  by  placing  a  series  of  obstructions  (Slrchelkin  spiral)  at  the  beginning  of  the  tube.  Obstructions 
(baffles)  in  GETF  experiments  played  a  similar  role  by  promoting  the  flame  acceleration  to  several  hundred 
meters  per  second  at  the  end  of  the  obstructed  section  of  the  tube.  For  all  three  cases  considered  here, 
the  flame  development  at  the  beginning  of  a  smooth  section  starts  from  a  fast  turbulent  flame  jet.  The 
subsequent  flame  development  eventually  leads  to  DDT,  which  occurs  many  tube  diameters  away  from  the 
beginning  of  the  smooth  section.  The  distances  Tddt  measured  from  the  beginning  of  the  smooth  section 
are  shown  by  open  points  in  Figure  5.1  and  approximated  by  line  1  as  a  function  of  tube  diameter  d. 

For  obstructed  tubes  (line  2),  the  initial  flame  is  ignited  by  a  spark  near  the  closed  end  of  the  tube,  and 
propagates  through  a  series  of  obstructions  until  DDT  occurs.  So  far,  there  are  only  two  experimental  data 
points  obtained  by  Kuznetsov  [72]  for  the  blockage  ratio  BR=0.3  and  the  distance  between  obstacles  S  =  d. 
Our  numerical  simulations  for  the  same  parameters  and  d  =  1.05  m  give  Lj^DT/d  =  17. 

Figure  5.1  shows  that  the  presence  of  multiple  obstructions  may  reduce  Ljjjjt  by  several  times.  A  rough 
extrapolation  to  d  =  3  m  (which  approximately  corresponds  to  the  mining  tunnel  size)  gives  L DDr/d  ~  30 
for  smooth  tunnels  with  violent  initiation,  and  ~  10  for  obstructed  tunnels  with  weak  spark  ignition.  In 
realistic  mine  conditions,  a  violent  initiation  may  occur,  for  example,  if  the  flame  is  ignited  and  accelerated 
in  a  heavily  obstructed  area,  and  then  enters  a  relatively  smooth  tunnel. 

These  estimations  are  valid  for  simple  geometries  and  uniform  obstacle  spacing  at  S  =  d.  Changing 
S  may  reduce  Anur/d-  For  example,  the  computations  performed  for  BR  =  0.3  and  d  =  105  cm  give 
LDDr/d  =  20  for  S  =  d  and  Loor/d  =  13  for  S  =  d/2.  Complex  geometries,  such  as  crosscut  tunnels, 
rough  walls,  mining  equipment,  or  partially  collapsed  areas  can  further  reduce  Tddt- 

Now  we  return  to  answer  the  overriding  question  of  whether  DDT  can  arise  in  coal-mine  conditions.  We 
have  shown  both  computationally  and  experimentally  that  it  is  possible  for  a  detonation  to  develop  from  a 
weak  spark  for  certain  NG-air  mixtures  in  certain  experimental  conditions.  Theoretically  and  computation¬ 
ally,  it  appears  that  DDT  will  develop  in  obstructed  and  smooth  tubes  over  a  range  of  NG  air  compositions. 
Preliminary  data  from  GETF  experiments  verify  that  DDT  is  achievable  in  obstructed  tubes  with  a  diameter 
of  1  m.  By  all  indicators  in  Figure  5.1,  these  results  will  extend  to  the  coal  mine  scale  where  tunnels  measure 
up  to  3  m  high,  6  m  wide,  and  are  hundreds  of  meters  long. 
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Figure  5.1  Run-up  distance  to  DDT  Ijddt  as  a  function  tube  diameter  measured  or  computed  by  different 
research  groups  in  smooth  (1)  or  obstructed  (2)  tubes  of  different  diameters  d.  GETF  data  (open  circles) 
are  obtained  in  experiments  with  9,  10,  or  11  baffles  placed  near  the  closed  end  of  the  tube  and  uniformly 
spaced  at  S  =  d  =  1.05  m.  Lddt  is  measured  from  the  last  baffle.  Cross  points  show  a  rough  extrapolation 
to  the  scale  of  a  typical  coal  mine  tunnel. 
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6.  Future  Work 


This  report  describes  the  NRL  research  efforts  to  answer  the  fundamental  safety  question:  Given  a  large 
enough  volume  of  a  flammable  mixture  of  NG  and  air,  such  as  may  exist  in  a  coal  mine,  can  a  weak  spark 
ignition  develop  into  a  detonation?  The  completion  of  the  four  specific  tasks  performed  to  address  this 
question  has  made  substantial  progress  towards  providing  an  answer. 

1.  Computations  that  solved  the  reactive  Navier-Stokes  equations  applied  to  DDT  in  stoichiometric 
methane-air  mixtures  showed  that  it  is  possible  to  calculate  the  complex  fluid  dynamics  of  reactive 
flows  in  obstructed  channels  and  simulate,  from  first  principles,  the  development  of  DDT.  Many  of 
the  physical  events  that  lead  to  DDT  can  now  be  calculated  and  are  well  understood.  Furthermore, 
the  simulations  agree  well  with  experiments  both  in  the  published  literature  and  recently  available 
with  the  GETF. 

2.  The  mechanics  of  DDT  in  methane-air  mixtures  contain  a  probabilistic  element  connected  with  the 
complexity  and  inherent  randomness  of  fluid  turbulence,  the  appearance  of  hot  spots  from  which 
detonations  develop,  and  other  features  of  combustion  processes.  The  computations  can  replicate 
the  stochasticity  of  the  DDT  process  and  the  variability  within  experimental  observations  of  DDT. 
(See  Section  3.2.2c  of  this  report.) 

3.  Combustion  processes  leading  to  DDT  in  fuel-lean  and  fuel-rich  methane-air  mixtures  can  be  com¬ 
puted  by  solving  the  reactive  Navier-Stokes  equation.  The  results  of  these  computations  agree  well 
with  experimental  observations  indicating  that  we  have  an  excellent  fundamental  understanding  of 
the  processes. 

4.  Computational  models  have  been  developed  and  tested  for  simplified  detonation  problems  for 
methane-air  mixtures  with  concentration  gradients. 

Recent  experimental  work  by  NIOSH,  performed  in  conjunction  with  theoretical  and  computational  work 
performed  by  NRL  has  extended  the  known  detonability  limits  for  methane-air  mixtures  to  the  5.3  to  15.5% 
range.  DDT  has  been  achieved  in  the  105  cm  diameter  tube. 

The  work  completed  so  far  by  NRL,  which  was  performed  in  conjunction  with  recent  experimental  work 
by  NIOSH,  has  made  substantial  progress  toward  answering  the  fundamental  underpinning  question  about 
DDT  stated  above.  We  have  not  yet  proved  conclusively  that  substantial  accumulations  of  NG  and  air  in 
coal  mine  tunnels  could  develop  into  detonation  when  ignited  by  a  weak  spark.  We  are,  however,  very  close 
to  providing  a  definitive,  scientifically  sound  answer  to  this  vexing  question. 

In  order  to  be  able  to  make  the  definitive  statement  on  DDT  in  coal  mines,  we  need  to  solve  three 
critical  problems.  We  have  to  be  able  to: 

(1)  Ensure  that  numerical  models  calibrated  on  small-scale  experiments  with  obstructed  channels  will 
give  the  correct  scaling  when  used  for  larger  geometries  typical  of  coal  mine  tunnels.  This  will  require  the 
model  validation  using  GETF  data  (GETF  experiments  have  not  yet  shown  DDT  in  the  obstructed  section 
of  the  tube),  and  possibly  performing  large-scale  experiments  in  LLEM  tunnels. 

(2)  Understand  and  predict  the  flame  behavior,  DDT,  and  detonation  propagation  in  mixtures  with 
spatially  nonuniform  compositions  that  are  likely  to  form  in  real-life  configurations.  Such  nonuniform  com- 
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positions  may  form,  for  example,  when  layers  of  methane  accumulate  at  the  top  of  a  horizontal  tunnel  without 
completely  mixing  with  air.  Solving  this  problem  is  within  capabilities  of  the  current  numerical  models. 

(3)  Understand  the  mechanism  of  DDT  in  a  smooth  section  of  a  channel.  This  is  one  of  the  key  unresolved 
problems  in  combustion.  This  type  of  DDT  is  observed  in  GETF  experiments,  but  is  not  reproduced  by 
the  current  models.  Addressing  this  question  will  require  a  development  of  more  sophisticated  numerical 
models,  possibly  with  a  detained  resolution  of  boundary  layers,  and  a  detailed  analysis  of  GETF  results, 
possibly  involving  additions  diagnostics. 

An  extended  NRL  and  NIOSH  effort  would  aim  to  produce  the  definitive  peer-reviewed  technical  papers 
on  metliane-air  detonation  and  DDT  applied  to  coal  mines. 
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